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SYMBOL INDEX 



Symbol Description 



a 


characteristic length, cylinder radius 


a 


dimensionless wavelength parameter, a = 27ra/L 


b 


dimensionless constant, H = eb 


c i 


th 

dimensionless wave force coefficient in the i 
direction 


d 


dimensionless relative cylinder depth 


d 


cylinder depth 


dC^ 


differential dimensionless arc length on the 
cylinder surface 


dC 2 


differential dimensionless length on the free 
surface 


f l 


first-order source strength function 


f 2 


second-order source strength function 


f* 


free surface particular solution portion of the 
second-order problem source strength function 


F( ) 


function of the quantities in parentheses 


F li 


dimensionless first-order wave force coefficient 
in the i^ direction 


F 2i 


dimensionless second-order periodic wave force 

th 

coefficient in the i direction 


F SS 
r 2i 


dimensionless second-order steady-state wave 

th 

force coefficient in the i direction 


g 


acceleration of gravity 


G 


Green's function 


G* 


modified Green's function for the particular 
solution portion of the second-order problem 
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Symbol Description 



h 


dimensionless relative mean water depth 


h 


mean water depth 


h. 

1 


first-order non-homogeneous boundary condition 

th 

function at the i nodal point 


H 


dimensionless relative wave height, H = H/2a 


H 


elevation difference between the wave crest 
and trough 


i 


complex plane portion, i = (-1) 2 


k 


wave number, k = 2 tt/L 


k i 


second-order non-homogeneous boundary condition 
function at the i nodal point 


K 


dimensionless Bernoulli constant 


K 


Bernoulli constant 


K 2 


second-order dimensionless Bernoulli constant 


K( ) 


function of the quantities in parentheses 


L 


dimensionless wave length 


L 


wave length 


m 


number of cylinder surface increments and 
nodal points 


n 


number of free surface increments and nodal 
points 


A 

n 


unit normal vector on the cylinder surface in 
the outward direction 


n ,n , 
x y 


spacial component unit normals in the horizontal 
and vertical directions, respectively 


0 


order of 


P 


dimensionless pressure coefficient 


P 


Pressure 


PV 


principal value 
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Symbol 



q 

r 
r ' 



t 



t 



U 1 

u 2 



u 



2 



U 



x,y 



x,y 



a 

6 

A 

V 

6 

£ 






n 



Description 
fluid velocity vector 
dimensionless linear distance 
dimensionless image linear distance 
real part of a complex number 
cylinder surface 
free surface 

dimensionless time, t = ct 
time 

first-order complex potential function 

second-order periodic complex potential 
function 

second-order time independent complex 
potential function 

non-periodic function 

dimensionless spacial variables in the 
horizontal and vertical directions, respectively 

spacial variables in the horizontal and 
vertical directions, respectively 

complex matrix 

complex matrix 

increment 

vector operator 

phase shift angle 

perturbation parameter 

dimensionless spacial variables corresponding 
in direction to x and y, respectively 

dimensionless free surface elevation 
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Symbol 


Description 


r\ 


free surface elevation 


e 


plane polar angle 


A 


half surface interval 


V 


dummy of integration variable 


V k 


positive roots of p^tan(iJ^h) - v = 0 


V 


dimensionless wave frequency 


7T 


constant, 3.14159 


p 


fluid density 


£ 


summation of terms 


0 


wave frequency 


4> 


dimensionless velocity potential 


$ 


velocity potential 


Subscripts 

i 


direction of a force or pressure component 
or nodal point location, depending on context 


j 


nodal point location 


n 


normal partial derivative 


t 


time derivative 


X, X 


spacial partial derivative in the horizontal 
direction 


y#y 


spacial partial derivative in the vertical 
direction 


1 


first-order component or x spacial direction, 
depending on context 


2 


second-order component or y spacial direction, 
depending on context 
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Symbol 



Description 



Superscripts 

I incident wave component 



o 



S 



homogeneous solution portion of second- 
order scattering problem 

scattering wave component 



particular solution portion of second- 
order scattering problem 
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I. INTRODUCTION 



The solution to the two-dimensional linear wave/structure 
interaction problem has been well-studied for both the finite 
and infinite depth cases. The fluid motion resulting from 
the interaction of a train of regular waves with a submerged 
horizontal circular cylinder in infinite depth water was 
first studied by Dean [Ref. 1]. Ursell [Ref. 7] later studied 
the problem anew, placing the solution on a rigorous basis. 

More recently, Ogilvie [Ref. 6] reconsidered the problem. He 
computed the first-order oscillatory forces and second-order 
steady-state forces for the following cases: (a) the cylinder 

held fixed, (b) the cylinder forced to oscillate sinusoidally 
in still water, and (c) the neutrally buoyant cylinder allowed 
to respond to wave motion. 

It can easily be shown that the steady-state part of the 
second-order forces can be computed from the first-order 
potential alone. Accordingly, Ogilvie did not obtain a 
complete second-order solution; the steady-state forces 
arise from the first-order velocity- squared terms in Ber- 
noulli's equation. In addition to the steady-state forces, 
second-order oscillatory forces also exist which have a 
frequency twice that of the first-order forces and represent 
the subject of the present work. 

This report reconsiders again the submerged horizontal 
cylinder problem, extending the analysis to include water 
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of finite depth. The cylinder is considered to be completely 
submerged and held fixed in a train of regular waves. The 
objective is to determine the complete second-order solution 
and, accordingly, determine the oscillatory second-order 
forces as well as the steady-state second-order forces. 

The problem is treated as a regular perturbation problem 
in the small parameter, incident wave height/ cylinder radius. 
Proceeding in this way the incident wave appears as a second- 
order Stoke 's wave. The boundary-value problems for both 
the first-order and second-order potentials are established 
and the solutions to both are obtained by use of the Green's 
function method. 
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II. THEORETICAL DEVELOPMENT 



A. FORMULATION OF THE PROBLEM 

A systematic representation of the problem considered is 
illustrated in Fig. (1) . A rigid right cylinder of radius a, 
submerged to a depth d in water of depth h, is fixed in 
a train of regular waves propagating in the positive x direc- 
tion. The basic problem is that of calculating the hydro- 
dynamic pressure on the immersed cylinder and the resulting 
force correct to the second-order in the wave height of the 
incident wave. 

Assuming the fluid to be irrotational , a velocity 
potential, §, may be defined as: 

q = V$(x,y,t) (1) 



/\ 

where q denotes the fluid velocity vector. (The barred 
quantities denote dimensional quantities.) Moreover, assuming 
the fluid to be incompressible, and homogeneous, the velocity 
potential must satisfy the Laplace equation, 

V 2 $(x,y,t) = 0 (2) 

within the fluid region. 

In addition to the differential equation, Eq. (2) , $ 
must satisfy certain boundary conditions. In specific, 
these are: 
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Figure 1. DEFINITION SKETCH 



1 . 



The velocity of the fluid normal to the bottom 



boundary must vanish. 



2. The velocity normal to the surface of the cylinder. 



defined by S^(x,y)=0, must vanish. 

3. The kinematic boundary condition on the free surface. 



defined by S 2 (x,y,t)=n (x,t)=0, must be satisfied. 
4. The dynamic boundary condition on the free surface 



must be satisfied. 



These boundary conditions may be stated mathematically as follows 



where n denotes the elevation of the free surface, g denotes the 
acceleration of gravity, and K denotes the Bernoulli constant. 

For convenience in carrying out the solution, the variables 
are next made dimensionless using the cylinder radius, and the 
wave frequency, a , as follows: 



<J>- (x , -h, t) =0 



(3) 



v$ • vs^o 



(4) 




(5) 



$-(x,n,t)+ [$-(x,n,t) + $-(x,n,t) ]+gn(x,t)=gK (6) 
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x/a 


y = y/a 


d = d/a 


h = h/a 


H/2a 


K = K/a 


2_ 

v = o a/g 


t = at 


a$/ga 


n = n/a 







( 7 ) 



H denotes the dimensionless wave height, K denotes the 
dimensionless Bernoulli constant, and t denotes the 
dimensionless time. 

Utilizing the parameters defined in Eq. (7) , the boundary 
value problem defined in Eqs. (2-6) may be rewritten concisely 
in dimensionless form as: 

7 2 $(x,y,t) =0 in the fluid region (8) 



<j>y(x,-h,t) = 0 (9) 

4> n (*,y,t) =0 on S 1 (x,y) = 0 (10) 

<J> x (x,n ,t) n x (x,t) - 4> (x,n,t) + vn t (x,t) = 0 (11) 

<t> t (x,n,t) + jU x (x,n,t) 2 + <}>y (x, n / 1) 2 ] + n(x,t) = k (12) 



B. PERTURBATION EXPANSION 

Expanding <}> and r| in terms of a perturbation parameter, 
e , provides a means for converting the nonlinear boundary-value 
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problem into a series of linear problems. The solution to 
each linear problem may be tractable whereas the nonlinear 
problem is unsolvable. 

Since <J>, n, and K are functions of the small parameter, 
e, they may be written in the power series: 



<j>(x,y,t;e) = X (x,y , t) (13) 

n=l n 



n(x,t;e) = X e n n (x,t) (14) 

n=l n 



and 



:=Xe r 

n=2 



K 



n 



(15) 



In Eqs. (11) and (12), <p contains n> and, therefore, e., 
implicitly. This can be converted to an explicit form in 
e by use of the Taylor series expansion: 



<p (x, n , t) 



00 , . s m 

J2 n ( x , t ? e ) 

m=0 mi 



r 9 m (j) (x,y , t) , 

L „ m J y=0 

3y * 



(16) 



The perturbation parameter, £, is related to the wave height 
in an, as yet, unknown manner. 

Equations (13-16) as well as appropriate derivatives 
similar to Eqs. (13-16) may now be substituted into the 
boundary-value problem given in Eqs. (8-12) to obtain a 
series of linear boundary- value problems for <f>^, etc. 
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That is, upon substitution and equating coefficients of like 
powers of e, the following problems for the first two terms 
in the expansion for c(» can be precipitated: 

First-order boundary-value problem (e): 



V 2 cf> 1 (x,y , t) = 0 




(17) 


$^y(x,— h,t) — 0 




(18) 


♦ ln (x,y,t) =0 on 


s 1 (x,y) = 0 


(19) 


4»iy (x , o , t) - vn lt (x,t) 


= 0 


(20) 


n 1 (x,t) + ^]_ t (x,o , t) = 


0 


(21) 



2 

Second-order boundary-value problem (e ) : 

V 2 <f> 2 (x ,y , t) = 0 (22) 

<fr 2 (x,-h,t) = 0 (23) 

^2n^ x,y,t ^ = 0 on s i^ x,y ^ = 0 (24) 

$ 2 (x,0,t) - vn 2t (x,t) = 

(25) 

-r^ (x,t) <t> ^ y (x,o , t) t <{>^ x (x, o , t) n^ x (x,t) 
n 2 (x,t) + < f> 2t (x,o,t) = -n-]_ (x,t) 4> lyt (x,o,t) - 
n 1 1 (x,t) <{) ly (x,0,t) - ^j(<t» lx (x,0 , t) 2 + <f> ly (X/O ,t) 2 ] + K 2 (26) 
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There is a striking similarity between the first-order 
and second-order problems; only the right hand side of the 
free surface boundary conditions differ. In the first-order 
problem the free surface boundary conditions are homogeneous 
whereas the second-order free surface boundary conditions 
are non-homogeneous , the right hand side being functions of 
the first-order velocity potential, first-order free surface 
elevation, and their derivatives. 

The first-order potential may be represented by a function 
which is periodic with the fundamental frequency. Accordingly, 
the complex potential, u^(x,y), may be expressed as: 

<j> 1 (x,y,t) = abR e [iu x (x,y) e -it ] (27) 

where R g denote the real part, a = 2iTa/L,L denoting the wave 
length, and b is an unknown real constant. 

C. THE FIRST-ORDER BOUNDARY VALUE PROBLEM 

Since the first-order boundary-value problem is linear, the 
total first-order potential may be expressed as the sum, 

= <!>i + 4>i (28) 

I S 

where <j>^ denotes the incident wave potential and <j»^ denotes 

the scatter potential which is due to the presence of the cylinde 

The space dependent part of the complex potentials are expressed, 

in view of Eqs . (27) and (28) as. 
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( 29 ) 



I ^ S 

+ U-^ 



The potential, u^, represents only a regular wave propa- 
gating along the x-axis and this potential must satisfy the 
first-order boundary-value problem when no rigid body is 
present. Therefore, <j>^ satisfying Eqs. (17-18) and (20-21), 
is given by 



I _ _ 1 cosh [a(y+h)] iax 
U 1 a cosh (ah) e 



(30) 



The parameters a and v are related through the well-known equation 
from linear wave theory. 

2 - 

v = -2—^ = a tanh (ah) (31) 

g 

which, in dimensional form, using a = ka, and k = 2tt/L is: 

a^= gk tank(kh) (32) 

Since the solution for u^ is known, the boundary-value 

s 

problem for u-^ may now be established. Substituting Eqs. (30), 

(29) and (27) into the first-order problem given by Eqs. (17-21) , 
and eliminating between Eqs. (20) and (21) results in a 
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boundary-value problem for the first-order scatter potential 
as f ollows : 

V 2 u^(x,y) = 0 (33) 

u^y(x,-h) = 0 (34) 

u ln (x ' y) = cosh 1 ( ah ' ) '[ n y sinh t a <y +h >) + in x cosh [a (y+h)] e lax 

on S 1 (x,y) = 0 (35) 

u^y(x,0) - vu^(x,0) = 0 (36) 

where n x and n denote the x- and y-components , respectively, 

/N 'N /V 

of the unit normal vector, n = in + jn , directed outward 

x y 

into the fluid. 

D. SECOND-ORDER BOUNDARY-VALUE PROBLEM 

In view of the linearity of the second-order problem, 
the same procedure as used in the first-order problem may 
be applied and leads to the representation of the second-order 
potential as the sum: 

4> 2 = < f >2 + 4>2 (37) 

where tj)^ denotes the second-order incident wave potential 

5 

and <J >2 denotes the scatter potential. 



24 



Considering the case where there is no body present, 
there is no scattered wave and therefore, <J>^ = ^ = 0 . 
Additionally, the boundary conditions on the immersed cylinder, 
Eqs. (19) and (24), are not applicable. Thus, when the substitu- 
tions of <J>^ for <j>^ and <j> 2 for <j> 2 are made into Eqs. (22) , 

(23), (25) and (26), along with Eqs . (20), (21), (27) and 

(31), and upon eliminating and n 0 between Eqs. (25) and 
(26) , the resulting boundary-value problem for the second- 
order incident wave potential becomes: 

V 2 4>2 (x,y,t) = 0 (38) 

4> 2y (x,-h,t) = 0 (39) 

<J>2 y (x,0,t) + v$ 2tt (x,0,t) = - |b 2 (a 2 -v 2 )R e [ie i2(ax-t) ] 

(40) 

The periodic solution to the incident^wave boundary-value 
problem specified by Eqs. (38-40) is known and may be expressed 
as 



4>2 = ^ab 2 vR e [iu 2 (x,y) e~ i2t ] 



(41) 



where the second-order complex potential, u 2 , is given by: 



u 2 (x,y) = 



1 cosh [2a (y+h) ] i2ax 

2a 1 ~ 4 . , . 
sinh (ah) 



(42) 
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Eqs. (41-42) are familiar in wave theory and are exactly the 
same form as that given by Ref. 4 for second-order Stokes' 
waves . 



viously determined solutions for the incident wave potentials, 
Eqs. (30) and (42), are utilized. After applying Eqs. (20), (21), 
(27), (29) and (41) and eliminating between Eqs. (25) and (26) 

the resulting boundary-value problem for the second-order scatter 
potential becomes : 



Having developed a solution for the second-order incident 

wave potential, <$>^ , the problem for second-order scatter 

potential may be determined. To this end,4>^=4>^ + <f>^ and 
I S 

<j >2 = <f >2 + 4>2 are substituted into Eqs. (22-26), and the pre- 



V 2 <j >2 (x,y, t) = 0 



(43) 




(44) 






4 

smh (ah) 




(45) 



on S 1 (x,y) = 0 





+ U(x) 



on y = 0 
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where U (x) is a non-periodic function generated by the 



substitution of <j>^ into Egs . (2.5) and (26). For brevity, 

this function is not written out since it is not needed in 

determining the second-order pressures and forces. 

g 

The boundary-value problem in <f > ^ given in Eqs , (43-46) 

is time dependent according to e l2t , except for the U (x) 

C 

term m Eq. (46) . Therefore, the solution for <j>. may be 
taken in the form 



$2 = f at>2vR e [ i [ u 2 ( x 'Y) e l2t + 1*2 (x,y)]] (47) 

where the last term denotes the time-independent portion of 

the complex potential. Separate boundary- value problems for 

S ~ s 

U 2 and U 2 arise from the substitution of Eq, (47) into Eqs. 

g 

(43-46) . However, as will be demonstrated only the <f> 2t term 

is required for determination of the pressure to the second 

order. Therefore, the time-derivative of the time-independent 

g 

part of <j >2 will be zero, and, accordingly, will not be developed 

g 

further here; the solution for U 2 only will be pursued. 

When Eq. (47) is substituted into Eqs. (43-46), the result- 
ing boundary-value problem for the second-order scatter potential 
takes the form: 

V 2 U2(x,y) = 0 (48) 



u 



S 

2y 



(x,-h)= 0 



(49) 
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u? (x,y) = \n sinh [2a (y+h) ] 

z sinh (ah) *- y 

+ in x cosh [2a (y+h) ]J e^ 2ax on S^(x,y)=0 
u 2y ( X/ 0 ) " 4 vu^ (x , 0) = f* (x) 



(50) 

(51) 



where 



, 2a r l SS IS , I S . IS I 

f*(x) = 37 t-u ly u lyy + u i u i y y - 6u ly u ly + v u ly u lyy 



(52) 



-IS - S 
- 4u lx u lx " 2u lx 



- 3u, ] n 

ly y=0 



There is considerable similarity in the first-order 
and second-order problems, the only difference in form being 
that the second-order free surface boundary condition, 

Eq. (51), is non-homogeneous . 

Because of the non-homogenity of Eq. (51), further use is 

S 

made of linear superposition defining U 2 as the sum. 



S S° , S* 

U 2 = ^2 + U 2 (53) 

S° • S * 

where U 2 denotes the homogeneous solution and U 2 the 

particular solution to the boundary-value problem as stated 

S * S° 

in Eqs . (48-51). More precisely, u 2 and u 2 are defined as 

the solutions to the boundary-value problems obtained from 
the substitution of Eq. (53) into Eqs. (48-51) as follows: 
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s * 

Particular Solution (u 2 ) : 

9 Q * 

V u* (x ,y) = 0 
u S 2y (x ,-h) = 0 



u 2y (x,0) - 4vu^ (x , 0 ) = f*(x) 



Homogeneous Solution (u^ ) 



2 s 

V u 2 (x,y ) = 0 



u 2y (x,-h) = 0 



s° s° 

u 2y (x,0) - 4 vu 2 (x , 0 ) = 0 



u~ (x,y) = j |~n sinh [2a (y+h) ] 

sinh (ah) L y 



+ in^ cosh [2a (y+h) ]Je^ 



2ax S , . 

- u 2n (x,y) 



on S^(x,y) = 0 



(54) 

(55) 

(56) 

(57) 

(58) 

(59) 

(60) 



By the division of the second-order problem into homogeneous 

S* 

and particular solutions, the non-homogeneous problem for u 2 
contains no boundary condition on the cylinder surface. The result- 

s * 

ing boundary-value problem for u 2 , as stated in Eqs. (54-56) , is 
identical to that associated with the linear problem for the potentia 
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resulting from a harmonic pressure variation of amplitude 
distribution f*(x) on the free surface in water of depth h. 

S° 

The homogeneous boundary- value problem for u^ , defined by 
Eqs. (57-60) , is similar in form to the first-order problem 
given by Eqs. (33-36) , the significant difference being the 
term 4v in place of v which occurs in the free surface boun- 
darv condition. Thus, the method of solution for U£ will 
be similar to that of the first-order scatter potential 

E. DEFINITION OF THE INCIDENT WAVE 

Having developed the first-order and second-order boun- 
dary-value problems, it is now appropriate to completely 
define the incident wave height and in the process determine 
expressions for the unknown constants b and . Solving 
Eqs. (21) and (26) for r>^ and r\^, respectively, and then 
substituting the results into the expression for the free 
surface elevation as given by Eq. (14) yields: 

n(x,t) = - e<t> lt + e t-* 2 t + *it*lyt + *ltt*ly * 61) 

' l ( +lx + +ly> + k 2 ! + 0<s3) 

where the potentials are evaluated at y = 0. Eq. (61) 
expresses the free surface elevation in terms of the total 
potential, and therefore, includes the effects of both the 
incident and the scattered waves. Hence, as it is of inter- 
est to obtain an expression for the free surface elevation 
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of the incident wave, the scatter potentials are set to 
S S 

zero, i.e. <f>^ = cj^ = 0. Evaluating Eq. (61) using the known 
solutions for the incident wave potentials, Eqs. (30) and 
(42), along with Eqs. (27) and (41) , then yields: 



D I (x,t) •= ebR e (e l(ax-t) 



0 rv2 , 2 2 . 

] + s 2 [S. (v -a ) 



4v 



+ K, 



(62) 



ab cosh (ah) [2 + cosh ( 2 ah) ] „ f i 2 (ax-t ),1 

+ 3 R le J 

sinh J (ah) 



+ 0(e J ) 



where n (x,t) denotes the dimensionless surface elevation 
for the incident wave with no cylinder present. If the x-axis 
is placed at the mean water line then the constant term must 
vanish and, therefore. 



K _ = 



_ b (a‘ 



2 , 

- v ) 



4 v 



(63) 



The remaining terms in Eq. (62) are time -depen den t f periodic 
functions, the second-order term having a frequency of twice 
that of the first-order. 

Defining the dimensionless wave height given in Eq. (7) 
as the difference in elevation between the crest and trough 
of the incident wave, eb must represent the wave amplitude. 
The second-order term in Eq. (62) having frequency twice the 
fundamental frequency makes equal contributions to surface 
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elevation at both the crest and trough, so that it contributes 
nothing to the wave height. Thus, in terms of dimensionless 
wave ‘height, b is given by 



H = eb = H/2a 



(64) 



where H denotes the elevation difference between the wave 
crest and trough. 

Expressing the incident wave profile in terms of the 
dimensionless wave height yields: 



It may be noted that Eq. (65) agrees with the expression 
for the second-order Stokes' wave given, for example, by 
Ref. 4. 

F. PRESSURES AND FORCES 

The pressure may be formulated by use of Bernoulli's 
equation arranged as follows: 

2 2 

P(x,y,t) = -P$=- - ip[$- + $- ] - pgy + pgK (66) 

where P(x,y,t) denotes the dimensional pressure and p denotes 
the fluid density. By use of Eqs. (7), (13-16), (27), (41), 



n 1 (x, t) = H cos (ax-t) 



(65) 




[2 + cosh (ah) ] cos [2 (ax-t)] 
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(47), and (64), Eq. (66) may be reduced to dimensionless form, 
and carried to the second-order in wave height. The result is 



-it. 



P(x,y,t) = -y - HaR e [u ] _e ] 



(67) 



H 2 a 2 L r6v 2 2 2, -i2t, 

— [ R e [ — U 2 “ u lx - U ly e ] 



+ ' U lxl + ' u ly 



2 x , 2 1 

I + ^ - 1 

a 



In Eq. (67) p(x,y,t) denotes the dimensionless pressure 
coefficient and is defined as: 



p(x,y,t) = P(x ^ ,t} (68) 

pga 

The first term in Eq. (67) represents the hydrostatic 
pressure as y is the dimensionless depth beneath the mean 
free surface (y = 0) . The second and third terms are harmonic, 
the third term having twice the frequency of the second, and 
representing the harmonic second-order contribution. The 
remaining terms in Eq. (67) are independent of time and result 
in the time-average or steady-state force components. 

Expressing the components of the wave force vectors in 
terms of integrals of the pressure over the cylinder surface 
area, the dimensionless force coefficients may be written 
as 



c i (t) 



S 



J P(x,y ,t)n i 
1 



dS 1 , 



i 



1,2 



(69) 
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where the dimensionless force coefficients in the x and y 
directions are defined, respectively, as 



c 1 (t) 



F x (t) 

pga 3 



c 2 (t) 



F^t) 

pga 3 



(70) 



(71) 



where F^(t) and F^(t) denote the force components. Addi- 
tionally, dS^ denotes a dimensionless differential arc 
length along the circumference of the cylinder. 

Applying Eq. (67) to Eq. (69) yields: 



C ± (t) = - y yn i dS 1 - Ha lt ]n i dS 1 



(72) 



^ra Co r 6v 2 2 2, -i2t 

" H [ 4V / R e [ — u 2 " u lx " u ly )e 



2 2 2 

+ ^ u lx^ + l u iyl + " 1 ^ n i ds i + °( h3 ) 



The first term in Eq{72) results from the hydrostatic 
pressure on the cylinder so it represents simply the buoyancy 
force. Disregarding the hydrostatic pressure the hydrodynamic 
force may be written as 
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(73) 



C i (t) = H F ±i cos (6 li ~t ) + H 2 

+ H 2 [F 2i cos ( <S2i-2t) + F 2 ?] + 0 (H 3 ) 



where the first-order, second-order, and steady-state force 
coefficients and phase shift angles are defined by comparison 
of Eqs. (72) and (73) as follows: 



i<5 



F 



2i 



2i 




a 

S 1 






)n.dS 

l 



1 



F 



SS 

2i 






l)n i dS 1 



(74) 



(75) 



(76) 



The dimensionless force coefficients, F^ and F 2 ^ are real. 

The first term in Eq. (73) represents the first-order 
forces of the fundamental frequency, a. The periodic portion 
of the second term in Eq. (73) represents the second-order 
contribution to the force at twice the fundamental frequency. 
The last second-order term represents the steady-state or 
time- independent contribution, sometimes called drift-force. 

Evaluation of the force coefficients defined by Eqs. 
(74-76) is the main objective of this work. Therefore, 
it is necessary to evaluate the first-order and second-order 
complex potentials, u^ and u 2 , as well as derivatives of u^ 
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on the surface of the cylinder. Additionally, evaluation 
of u^ and its derivatives on the mean free surface (y = 0) 
is required in order to carry out the solutions. 
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III. METHOD OF SOLUTION AND NUMERICAL PROCEDURES 



A. SOLUTION OF THE FIRST-ORDER PROBLEM 

The use of a Green's function to express the solution 
to the first-order boundary-value problem was first formu- 
lated by John [Ref. 5] , and applied to submerged ellipsoids 
by Garrison and Rao [Ref. 3] . This method is considered to 
be applicable, in principle, to arbitrary shapes and is 
mathematically the most straightforward. The Green's func- 
tion, G, of unit strength which satisfies Eq. (33) as well 
as the boundary conditions, Eqs. (34) and (36) is given by: 



G (x,y; 5,n; v) 



In r 



i n t-i + op V f f cosh [ y (y+h) ] cosh [y (n+h) ] 

qJ L (coshyh) (vcoshyh-ysinh yhT 



-yh sinh(un) sinh(uy) 

sinE ( yh) 



cos I x-C j dy 



(77) 



- l 



4tt 

2a^h+sinh (2a^h) 



cosh [a^ (y+h) ] cosh [a^ (y+h) ] cos [a^ | x~C 



] 



where : 



r 



r ’ 



[ (x-C) 2 
[(x-C) 2 



+ (y-n ) 2 ] 2 



+ (y+y) 2 ] 2 



(78) 

(79) 



The symbol, a^, is defined in terms of h and v as the 
solution to the equation 
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F (a^,h,v) 



0 



(80) 



where F is defined by 

F(a^,h,v) = a^ tanh(a^h) - v (81) 

Comparison of Eqs . (31) , (80) and (81) demonstrates that a^ 

is clearly the equivalent of a. However, this will not be 
the case for the second-order problem, thus, requiring the use 
of separate notation. In Eq. (77), PV denotes the principal 
value of the integral. 

The Green's function given in Eq. (77) was also given in 
series form by John [Ref. 5] as 

2 2 

“ (y k + v )cos[y. (y+h) ] -y. |x-$| 

G(x,y;£,n?v) = 2 tt L, — * cos [y. (n+h) ] e 

k=l v l J k " h Vi k ("lj c + v ) 

ia il x -C| 

4ircosh [a 1 (y+h) ] cosh [a, (n+h) ] e 

- i (82) 

2a^h + sinh (2a^h) 

where a^ is as defined in Eq . (80) and the quantities y^ 

are defined as the real positive roots of the equation 

K(y k ,h,v) = 0 (83) 

where the function K is defined as 
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K(p,h,v) = y tan(yh) + v 



(84) 



s 

Followxng the Green's function method of solution, u^ 



is written as the integral over the cylinder arc length, 
S as: 



where (£,n) denotes points on the immersed surface, f^(?,n) 
denotes the unknown source strength, and dS^ = dS^/a denotes 
the differential arc length on the cylinder surface, made 
dimensionless with the characteristic length, a. 

Although the solution to the first-order boundary-value 
problem as stated in Eqs. (33-36) is given by Eq. (85) , the 
source strength function, f^(£,n)/ must be determined in 
order to evaluate the potential. From potential theory, 

S 

the normal derivative of the potential, u^(x,y), on the 
surface of the cylinder is 



where G (the normal derivative of G) may be determined by 
n 

differentiation of either Eq. (77) or (82) in a straight- 
forward manner. 

Applying the final boundary condition, Eq. (35) , leads 
to an integral equation which may be solved for f^. 




dS 



1 



(85) 




( 86 ) 
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2 f l (x,y) + 2 it 



(87) 




(S*n) 



G n (x,y; £,n,v) 



dS 



1 



cosh (ah) 



jn^sinh [a (y+h) ] + in^cosh [a (y+h) ]1 1 



lax 



The solution for from Eq . (87) may then be used in Eq . 

S 

(85) to determine the potential, u^ . 

B. SOLUTION OF THE SECOND-ORDER PROBLEM 

The similarity in form of the boundary-value problem 
for the homogeneous part of the second-order potential, 

Eqs. (57-60), to the first-order potential is now utilized. 
Since the only significant difference occurs in Eq. (59), 

S° 

where v is replaced by 4v, we may represent u 2 in a form 
similar to Eq. (85) as 

s ° if 

u 2 (x,y) = /f 2 (£j,n)G(x,y;£,v;4v) dS 1 (88) 

S 1 

where G(x,y; £,n; 4v) is defined by replacing v by 4v in Eqs. 
(77) and (82). Therefore, a 2 is defined by 



F (a 2 ,h , 4v ) = 0 (89) 

and a 1 is replaced by a 2 in Eqs. (77) and (82), also. In the 
case of Eq. (82) , is defined as the real positive roots of 
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K (y k ,h,4v) 



0 



(90) 



The source strength function, f 2 (£,n)/ appearing in 
Eq. (88) is again determined by applying the kinematic 
boundary condition on the cylinder surface as given by Eq . 

(59) . The integral equation for f 2 may then be given by 

\ f 2 (x,y) + ~ Jt 2 (£,n)G n (x,y;£,n;4v) ds 1 (91) 

S 1 

”n sinh [ 2a (y+h) ] + in cosh (2a (y+h) ] 

— y x 

g 

sinh (ah) 
on S (x,y) = 0 

S * 

However, to solve Eq. (91) for f 2 , first a solution for u 2 

S * 

must be obtained and u 2n evaluated on the cylinder surface. 

S * 

Thus, at this point the solution to u 2 is sought. 

The complex potential u 2 is defined as the solution to 
the boundary-value problem, Eqs . (54-56). The form of this 

problem is recognized as being equivalent to that of fluid 
motion produced by a periodic pressure distribution (as 
specified by f*(?)) with frequency 2a on the free surface. No 
cylinder is considered to exist in the fluid region. 

The solution to this problem may be formulated as an 
integral over the free surface extending from negative 
infinity to positive infinity. Using the present 



i2ax S* . 

e “ u 2n (x ,y ) 
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( 92 ) 



S* 



dimensionless representation, becomes: 

u 2 * ( x , y ) = ^ s ^f* (?)G* (x,y; C/0;4v) d£ 

2 

The Green's function, G* , for this problem is given in Ref. 8 
Upon transforming the solution to the non-dimensional form 



G*(x,y;5,0,4v) = -2PV / coshtu (y+h)]cos tii(x-e) Idu 

-°° 4vcosh(yh) - ysinh(yh) 



4ircosh ^h] cosh [a 2 (y+h) ] cos [a 2 (x-£) ] 



-l 



2a2h + sinh(2a2h) 



(93) 



where a 2 is defined by 



F(a 2 /h,4v) = 0 



(94) 



The normal derivative of is required in order to 

completely define the kinematic boundary condition on the 
cylinder as given in Eq. (60) . From Eq. (92), this becomes 



u 



S* 

2n 



(x,y) = i /f*(?)G* (x,y;C,0;4v) d£ 

C J 



(95) 



The derivative of G* in the direction normal to the cylinder 
surface may be obtained by differentiation of Eq . (93) in a 

straightforward manner. Thus, using the values for f*(£), 
as defined in terms of the first-order solution only, and as 
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S* c * 

given in Eqs. (52) and (56), ^ and U 2 n may be determined. 

S * S° 

Using U 2 n , U 2 may then be determined to complete the 
second-order boundary-value problem solution. 

C. NUMERICAL METHODS 

Determination of the forces on the cylinder surface 

requires solving for the potentials u^, U 2 and the derivatives 

of u-j^ at points on the surface, as shown in Eq. (72) . Thus, 

the problem is now one of solving for both the first-order 

and the second-order scatter potentials, and requires solving 

for the source strength functions f^ and as well as the 

free surface pressure distribution source strength function, 

* 

f . In view of the complexity of the equations it is natural 
to attempt a numerical solution. 

g 

The first-order scattering potential, u^ , as given by 

Eq. (85) is dependent upon the first-order source strength 

function, f ^ . Thus, it is necessary to solve Eq. (87) for 

g 

f^ to determine u^ . A numerical solution may be developed 
by dividing the cylinder surface into elements of length 
A6 = 2^/m, with the center of each element assigned an index. 
Since f^(£,ri) is recognized as a well-behaved function for a 
smooth surface cylinder, it is appropriate to define parameters 



ct. . (v) = i 

13 IT 



As 



y^n^ X i ,y i ; ^ ,n;v) dS l 



lj 



i , j 1,2,3, ...,m 



(96) 
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cosh (ah) 



1 



[ny(x i ,y i ) sinh [a (y^+h) ] 



(97) 



-■ iax. 

+ in (x.,y.) cosh [a (y . +h) ] e 

X 11 1 




(98) 



and, thus, Eq. (87) may be approximated by the complex matrix 
equation. 



cylinder . 

For the purpose of evaluating a and 8, either Eq . (77) 

or Eq . (82) may be used. For evaluation of a and 3 when i = j 

Eq. (77) must be used to allow separate treatment of the 
logarithmic singularity as r approaches 0. The difficulty 
encountered with the logarithmic singularity may be overcome 
by carrying out its integration analytically. Garrison 
[Ref. 2] showed for the calculation of that the contribu- 

tion of the normal derivative of the In (r) in Eq. (77) 
integrated over the singular element of A0 is A 0/2, not only 
for the nodal point (x^,y^) = (£,n)/ but for all nodal points - 




i/j = 1,2 



(99) 



Upon evaluation of cuj (v) using Eq . (96), the inversion of 

Eq. (99) may be carried out on the digital computer to 
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on t ^ le c Y linder surface. Additionally, to calculate 



integral of the In (r) over the singular element of A0. 

A second singularity arises in the evaluation of the 
infinite integral occurring in Eq. (77) . The first term in 
the integrand is singular at y = a. Since the numerator 
approaches zero like (y-a) near a. Garrison [Ref. 2] sub- 
tracted l/(y-a) from the singular term in the range 0 <_ y £ 2a 
and added the contribution of the integral of 1/ (y-a) (which 
in this case was zero ) . This converted the singular integrand 
to a regular function. Applying this technique to numeri- 
cally evaluate the resulting integral between 0 and 2a, 
Simpson's three-eights rule is used with an odd number of 
subdivisions, thus ensuring that the point y = a is not 
encountered . 

With f^ determined from the inversion of Eq. (99), the 
first-order potential and its derivatives may be evaluated 
on the cylinder surface. Replacing the surface integrals 
with summations, Eq. (85) and its derivatives may be written 
as 





i, j = 1,2,3 



( 100 ) 




( 101 ) 




( 102 ) 
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c c 

where , u iyi' etc • denote functions evaluated at the 
i th nodal point on the cylindrical surface. The complex 
matrices in Eqs . (100-102) are defined as follows: 



e ij( v) 



B xij (v) 



6 • • (v) 

yin 



_1_ 
2 TT 



AS 



yG(x i ,y i ; 5,n;v)dS 1 i,j = 1,2, ...,m (103) 



lj 



_1_ 

2 7T 



AS 



^G^ (x^ / ^ ' 0 ? v) dS^ 



(104) 



ID 



= /g (x . ,y . ; £ ,r) ; v) dS-, 

2TT A cyy i i 1 



(105) 



As lj 



The first-order incident potential and its derivatives may 
be evaluated directly at each nodal point after differentia- 
tion of Eq. (30) as needed, and combined with the scattering 
potential and its respective derivatives. 

To solve the second-order problem, the function f*(£) 
required in Eq. (95), and defined by Eq. (52), must be 
evaluated numerically . Dividing the mean free surface 
(y = 0) in the vicinity of the cylinder into n equal incre- 
ments, f* is evaluated at each nodal point. The numerical 
solution again uses Eqs. (100-105) with y^ = 0 for the 

5 

purpose of evaluating u 1 ^ and its derivatives at the mean 
free surface nodal points. The first-order wave potential 
and its derivatives are evaluated directly at each nodal 
point, and combined with the scatter potential and its 
derivatives to solve Eq. (52) for f*. 
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Having determined f* at all nodal points on the free 

S* 

surface, the boundary-value problem for U 2 given by Eqs. 

S* S* 

(54-56) may be solved, i.e. u,, and u 0 may be evaluated 
on the cylinder by solution of Eqs. (92) and (95) respec- 
tively. Expressing the integrals as numerical summations in 
complex matrix form yields: 

U 2 ni = f^at.(4v) i = 1,2, ...,m (106) 

J J j = 1,2, . . . ,n 

u 2i = f j 6 ij (4v) (107) 

where 



* 

a i j (4v) 
6* j (4v) 




Y i ; C/Ti; 4v) 



5/0 ; 4v) 



dC 



d£ 



(108) 



(109) 



s* s* 

Thus, using f * , both u 2 and u 2n may be directly evaluated 
at each cylinder surface nodal point. 

To solve the boundary-value problem for the second part 

s ^ 

of the second-order scatter potential, u 2 , as specified 
by Eqs. (57-60) , Eq. (88) must be evaluated. However, to 

S° 

evaluate Eq . (88) for ^ , the second-order source strength 

function, m ust be evaluated by numerically solving 
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the integral equation, Eq. (91) . Replacing the integral 
equation with a complex matrix summation: 

f 2 £ + f2j a j_j(4v) = 2k^ i,j - 1,2, ...,m (110) 

where a^j (4v) is defined by Eq. (96) , 



f 2 j = f 2 < x j'yj> 



and, 
k. = 



= i [n (x. ,y. )sinh[2a(y.+h) ] 

sinh (ah) *- y 



( 111 ) 



( 112 ) 



i2ax . 



1 g* 

+ in x (x i ,y i )cosh[2a(y i +h) ]J e - ^ 2 n (x i ,y i ) 



Once a -j_j (4v) is evaluated, Eq. (110) may be solved to obtain 
the second-order source strength function, f 2 , at each nodal 
point on the cylinder surface. Stating Eq. (88) in summation 
form, 




f 2j e ij 



(4v) 



i , j 1,2,3, ...,m 



(113) 



where S-^(4v) is defined by Eq. (103). Using Eq. (113), 

S° 

U 2 may be determined at each nodal point on the cylinder 

S* 

surface, then combined with and the second-order inci- 

dent potential, u* , to evaluate the total second-order 
potential, u 2 • 
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Since the first-order potential, u^, and its derivatives, 

and u^ , as we ll as the second-order potential, Uj / are 
now known at each nodal point on the surface, the first- 
order, second-order and steady-state force coefficients and 
phase shift angles may be determined using Eqs . (74-76). As 

shown earlier, the integrals may be replaced by summations, 
using nodal point values and the arc length increment, A6. 

Once each integral is evaluated, then the force coefficient 
becomes the absolute value or modulus of the complex integral 
result and the phase shift angle is the angle whose tangent 
is the imaginary part over the real part of the complex 
integral result. 

D. COMPUTER SOLUTION 

In order to utilize the method of solution described 
in Sections B and C, a computer program was developed to 
carry out the indicated calculations. Although accuracy was 
a prime consideration, an equally important requirement was 
to minimize computer run time. To achieve this, every attempt 
was made to utilize symmetry and to generate certain constants 
and matrices to eliminate redundant computations. 

The program consists of a main program that flows in the 
order of computation presented in Sections B and C, along 
with four subroutines. Two subroutines to solve the first- 
order Green's function, using both forms of the function, 

Eqs. (77) and (82) , were developed by Garrison [Ref. 2]. 

These subroutines, GREEN and GREENS respectively, have been 
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modified to include calculation of the first-order scatter 

potential derivatives, evaluation of the first-order scatter 

potential and its derivatives on the mean free surface, y = 0, 

and evaluation of the modified Green's function, G* , for 

S * S * 

the determination of ^ and U 2 n at each cylinder nodal 
point. 

For elements of the a and $ matrices corresponding to 
small values of the parameter (x-£) , GREEN is used while for 
larger values of (x-£) , GREENS is used. With the exception 
of the diagonal elements in Eqs. (96) and (103-105) and a few 
surface nodal points in Eqs. (108-109), the majority of 
elements of a and 6 are calculated by GREENS. This is most 
fortunate as the series form converges rapidly, requiring much 
less computer time than the equivalent integral form. 

Subroutine GEODAT reads the input geometrical data, 
generates the matrices h^ and as defined by Eqs. (97) and 
(112) respectively, and calculates certain geometrical 
parameters and matrices for repeated use in GREEN and GREENS. 
Subroutine COMAT inverts complex matrix equations and, there- 
fore is used to invert Eqs. (99) and (110) to determine f 1 and 
respectively. 

A cross-reference between text and computer program 
nomenclature is given in Table I . 
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TABLE I: Computer Program — Text Symbol Cross-Reference 



Text 


Computer 

Program 


Text 


Computer 

Program 


a 


A 


n 


ANX ( I ) 






X 


d 


D 


n y 


ANY (I) 


f l 


F (1 , 1) 


u i 


Ul (I) 


f 2 


Fl (1,1) 


U lx 


U1X (I) 


f* 


FS (L) 


u iy 


UIY (I) 


i — ! 
i — 1 


Cl(l) 


u* (surface) 


U1IS 


F 12 


Cl (2 ) 


u lx ^ sur ^ ace ) 


U1ISX 


F 21 


C2 (1) 


u F y (surface) 


UlISY 


F 22 


C2 (2) 


u^^ (surface) 


UlISYY 


F SS 

21 


C 3 ( 1 ) 


s 

u^ (surface) 


U1SS 


F SS 
r 22 


C3 (2) 


s 

u lx ^ sur ^ ace ) 


U1SSX 


G 


GI J , GI JEXT 


s 

u^ (surface) 


U1SSY 


G* 


GI J , GI JEXT 


s 

u iyy ( surface ) 


U1SSYY 
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IV. DISCUSSION AND RESULTS 



A. SELECTION OF COMPUTER PROGRAM PARAMETERS 

The computer program was developed on the IBM-360 
computer using FORTRAN H. All subroutines were compiled on 
DATA CELL to minimize compile time and core size. In order 
to maximize accuracy while minimizing computational time, 
two key computer parameters were evaluated to determine 
optimum values: cylinder nodal points and free surface model 
points . 

Using the first-order solution only, the number of nodal 
points on the cylinder surface was varied from 12 to 60 and 
compared with the results determined by Garrison [Ref. 2], 
The minimum number of cylinder surface nodal points that 
provided accuracy within one percent of those obtained by 
Garrison using 60 points, was 24. Accordingly, the good 
correlation from 60 points down to 24 points led to the 
selection of 24 nodal points on the cylinder surface. 

Since the free surface integral was to be carried out 
from to + 00 , the next step was to establish the finite 
interval of convergence on the surface as well as the sub- 
division size. After the outer limits of the free surface 
integral were determined, the subdivision size was varied 
from 4 to 128 subdivisions per second-order wave length 
holding the total interval constant. Above 64 subdivisions 
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the results did not vary more than one percent, leading to 
the selection of 64 subdivisions per second-order wave length 
on the free surface for purposes of making calculations. 

Additionally, the total interval, 2A, was varied from 
one to eight wave lengths in both the positive and negative x 
directions. The second-order force coefficients and their 
respective phase shift angles converged to a small value 
varying periodically. Convergence to this limit cycle 
occurred between the first and second wave length from 
the center of the cylinder. Therefore, a total interval of 
from -2 to +2 second-order wave lengths was chosen as adequate 
for generation of the numerical results. 

Although the values of the second-order force coefficients 
and their respective phase shift angles tended to converge 
with increasing the limits of the free surface integration, 
there remained, in general, the small limit cycle as noted 
above. The amplitude of the cycle was a function of the 
parameter, a, and was of significant magnitude only for 
values of a less than 0.25. The effects of a on the limit 
cycle for one of the second-order parameters, horizontal force 
coefficient, is demonstrated in Fig. (2) . This represents a 
plot of percent variation of the force coefficient from the 
mean versus the surface interval size, A/(L/2), corresponding 
to values of a of 0.25, 0.20 and 0.15. Figure (2) is 
representative of the behavior of all second-order force 
coefficients and phase shift angles and demonstrates the 
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SURFACE INTERVAL \/(L/2) 

Figure 2. SECOND- ORDER HORIZONTAL WAVE FORCE 
COEFFICIENT VERSUS SURFACE INTERVAL 

d = 2.0, h =5.0 



increase in the amplitude of the variation with a decrease 
in a. The amplitude of the periodic variation ranges from 
up to fifty percent at a wave number of 0.15 and depth of 
d = 1.4, to less than two percent for all values of wave 
number above 0.25. 

B. RANGE OF APPLICABILITY 

In a complex computer program of the type developed to 
carry out the present numerical scheme, certain limits with 
respect to numerical stability, etc. show up. These computing 
limits arise for various reasons such as overflows caused by 
computing hyperbolic functions of very large numbers, 
numerical convergence instabilities, etc. While no attempt 
was made to investigate each of these numerical problems , a 
list of the limitations of the particular program developed 
in this work to carry out the computations are as follows: 

minimum a — dependent upon acceptable error, but 0.25 
or less 

maximum a — deep water condition reached 

minimum h — 3 

maximum h — deep water condition reached for a > 0.25, 
h = 20 for a < 0.25 

minimum d — from 1.4 at h equal to or less than 5 down 
to 2.5 at h = 20 

minimum SMIN — 0.12 when cylinder depth is other than 

near the free surface 

maximum SMIN — 0.30 when cylinder depth is near the 

free surface 
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C. RESULTS 



A representative set of results for the first-order and 
second-order force coefficients have been generated for a 
water depth of h = 5.0, submergence depths of d = 1.5, 2.0, 

2.5, and 3.0, with a ranging from 0.15 to 1.2. All of the 
numerical results presented herein were based on 24 sub- 
divisions on the circular cylinder and 64 subdivisions per 
second-order wave length on the free surface integration. 

The surface integral was carried out from x = -L to +L since 
this was found to be adequate for convergence. Since the 
second-order wave length is half the first-order wave length, 

L, the total interval, is equal to four second-order wave 
lengths . 

The first-order horizontal and vertical force coefficients 
are presented in Figs. (3) and (4), respectively. It may be 
noted that, in general, the forces decrease with depth of 
submergence according to expectation since the wave action 
dies out with depth. 

The second-order horizontal and vertical force coefficients 
are presented in Figs. (5) and (6), respectively. Generally, 
the results show the second-order effect to be relatively more 
important as the cylinder approaches the free surface. It 
might be suspected from this that the second-order contribution 
would be much more significant in the case of surface piercing 
or floating bodies. 
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The horizontal and vertical steady-state force coefficients 
are shown in Figs. (7) and (8) , respectively. These results 
show that the horizontal steady-state force coefficient is very 
small in general. The horizontal force can be shown to be 
proportional to the momentum flux of the reflected wave and 
the reflected wave is small. In fact, in the case of infinite 
water depth. Dean [Ref. 1] showed that the reflection was 
exactly zero, and, accordingly, the steady-state horizontal force 
is likewise zero. 

It may be noted that the steady-state vertical force is 
positive. This results from the fact that the velocities are 
largest on the top of the cylinder and, accordingly, the average 
pressure is reduced on the top of the cylinder. 

The phase shift angles of the first-order and second-order 
forces are shown in Figs. (9) - (12) . 

To demonstrate the effect that the second-order terms have 
on the forces on the cylinder and their respective phase shift 
angles, a comparison was made between the first- and second- 
order results. Figures (13) - (13) are plots of the horizontal 
and vertical dimensionless forces on the cylinder versus time 
over one complete wave cycle for both the first-order solution 
alone and the combined first-order and second-order effects. 
Additionally, a plot of the incident wave is included to 
demonstrate the phase shift involved (the amplitude of the 
incident wave has no significance) . A mean water depth, h, of 
5.0, a cylinder depth, d, of 1.5 and a wave height, H, of 
0.5 were used, with the wave number, a, assigned three values, 
0.25, 0.5 and 1.0. 
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0 0.2 0.4 0.6 0.8 1.0 1.2 

WAVE NUMBER, a = 2ir~a/L 

Figure 3. FIRST-ORDER HORIZONTAL WAVE FORCE COEFFICIENT • 
h =5.0 
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WAVE NUMBER, a = 2ira/L 

Figure 4. FIRST-ORDER VERTICAL WAVE FORCE COEFFICIENT : 
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0 0.2 0.4 0.6 0.8 1.0 1.2 

WAVE NUMBER, a = 27ra/L 

Figure 5. SECOND-ORDER HORIZONTAL WAVE FORCE COEFFICIENT ; 
h = 5.0 
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WAVE NUMBER, a = 2va/L 

Figure 6. SECOND-ORDER VERTICAL WAVE FORCE COEFFICIENT = 
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0 0.2 0.4 0.6 0.8 1.0 1.2 

WAVE NUMBER a = 2 tto/L 

Figure 7. STEADY- STATE HORIZONTAL WAVE FORCE COEFFICIENT: 
h =5.0 
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WAVE NUMBER, a = 2 tto/L 

Figure 9. FIRST-ORDER HORIZONTAL PHASE SHIFT ANGLE 
h = 5.0 
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0 0.2 0.4 0.6 0.8 1.0 1.2 

WAVE NUMBER, a =27ra/L 

Figure 10. FIRST-ORDER VERTICAL PHASE SHIFT ANGLE : 
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Figure 14. VERTICAL WAVE FORCE: a =0.25, h = 5.0, d = 1.5, H = 0.5 
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0 1.0 2.0 3.0 4.0 5.0 6.0 

TIME t (rad) 

Figure 15. HORIZONTAL WAVE FORCE: a = 0.5, h = 5.0, d = 1.5, H = 0.5 
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F inure IF \/FRTIFAI VA/AX/F FHRCF : n = r\R 
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Figure 17. HORIZONTAL WAVE FORCE= 



Second-Order 
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Figure 18. VERTICAL WAVE FORCE : 



V. CONCLUSIONS 



A computer program has been developed to carry out the 
numerical solution of the second-order wave — cylinder 
interaction problem. The first-order, second-order, and 
steady-state force coefficients were determined for the 
submerged horizontal cylinder. 

An approximate method for dealing with the second-order, 
nonhomogeneous free surface boundary condition was developed 
which appears to converge except at small values of a, i.e., 
at very large wave lengths. 
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APPENDIX A 

COMPUTER PROGRAM. LISTING 

THIS PROGRAM CALCULATES WAVE FORCES AND THEIR PHASE 
SHIFT ANGLES FOR WAVE INTERACTION WITH A SUBMERGED 
HORIZONTAL CYLINDER 

A = 2*PI*CYLINDER RADIUS/WAVE LENGTH = SIGMA**2* 
CYLINDER RADIUS/ACCELERATION OF GRAVITY 
H = MEAN WATER DEPTH/CYLINDER RADIUS 
D = CYLINDER DEP TH /C Y L I NDER RADIUS 

NPTS = NUMBER OF CYLINDER SURFACE ELEMENTS AND NODAL 
POINTS FOR NUMERICAL EVALUATION 

NSPTS = NUMBER OF FREE SURFACE ELEMENTS AND NODAL 
POINTS FOR NUMERICAL EVALUATION 



COMPLEX GI J » GN I J » GNJ I ,GX I J , GX J I , GY I J , GY J I , GI JEXT, GYY, 
1DET,SUM1,SUM2,SUM3,SUM4,U1I S , U II SX , U 1 1 S Y ,U 1 I SYY , 
2U1SS»U1SSX,U1SSY»U1SSYY 
COMPLEX ALPHA (24, 2 A) , BETA( 24,24) , BETAX( 24,24) , 

1BETAY ( 24, 24) ,HH( 24,1) ,F(24, 1) , PK(24, I) , FI (24,1) , 

2F S ( 500 ) , U 1 ( 24 ) » U 1 X ( 24 ) ,U1Y(24) ,U2(24) ,J2SC1(24) , 
3U2SCNI ( 24) ,U2SCO( 24) 

DIMENSION X<24) , Y( 24) ,ANX( 24) , A NY ( 24 ) , C HY ( 2 5 ) ,CHY2( 25) 
1, SHY (25) , SHY 2(25) ,C0EFG(200) ,COEFG2( 200) , 

2COSAMU ( 200 ,2 5) ,C0SAM2( 200,25) , SI NAMU( 2 JO ,25 ) , 

3 S I NAM2 ( 200,25) ,AMU(2C0) ,AMU4(200) ,SH2Y( 24) ,CH2Y( 24) 
DIMENSION C1(2),C2(2),C3(2), PHASE 1(2) , PHASE 2( 2) 

COMMON/ CP X /HH » P K 
CCMMON/VAR/X, Y, ANX, ANY 

C0MMCN/GSHY/ChY,CHY2 , SHY, SHY2 , SH2Y , CH2Y 
C CMMON/GMU/ COSAMU ,C0SAM2 , S I NAMU , S INAM2 , AMU , AMU 4 , COEFG, 
1C0EFG2 

CCMMON/CONST/H»D» DELTHE ,SMI N,NPTS , NSPTS 

COMMON/CPI/PI 

EQUIVALENCE (Hh(l),F(l)) 

EQUIVALENCE ( PK ( I ) , F 1 ( I ) ) 

PI=3. 141592 



READ INPUT DATA AND CALCULATE REQUIRED REPEATING 
DATA AND ARRAYS 



CALL GEODAT ( A , A2 , ANU, ANU4 , SH2 AH , SH2AH2 , SHA H , SHAH2 , 
1CHAH,CHAH2,A0,AA,3B,CC,DD,A02,AA2,BB2,CC2,DD2) 

DO 100 1 = 1, NPTS 
DO 100 J=1,I 
XV=X( I ) 

Y V=Y ( I ) 

XC=X( J ) 

Y C=Y ( J ) 

ANX 1 = ANX ( I ) 

ANY 1 = ANY ( I ) 

ANXJ=ANX( J ) 

ANY J=ANY( J ) 

I F ( ABS ( X ( I)-X(J) ) .LT.SMIN) GO TO 50 
SHY I = SH Y ( I ) 

CHY I = CHY ( I ) 

SHY J= SH Y ( J ) 

CHY J=CHY ( J ) 

CALL GREENS l A , ANU , S H2 AH , SH AH , C H AH, COS AMU , S IN AMU , AMU , 
1C OEFG , S HY I , CHY I , SHY J ,CHY J, I,J,XV,YV,XC,YC,ANXI,ANYI, 
2ANXJ»ANYJ»GIJ,GNIJ, GNJ I»GXIJ,GXJI,GYIJ,GYJI , GYY ) 

GC TO 90 
50 CONTINUE 

CALL GREEN ( A , ANU , SH2 AH, SHA H , CH AH , AO , AA , BB , CC , DD , I , J , 
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1XV,YV,XC, YC, ANXI , ANY I , GI J , GN I J , GX I J , GY I J , GY Y ) 

IF(I.EQ.J) GN J I =GNI J 
IF (I.EQ.J) GX J I =GX I J 
IF (I.EQ.J) GY J I = GY I J 
I F ( I.EQ.J) GO TO 90 

CALL GREEN < A , ANU , SH2 AH , SHA H, CHAH, AO , AA , 86 , CC , DD , J , I , 
1XC,YC,XV,YV,ANXJ , ANYJ,GI JEXT, GNJI ,GXJI , GYJI , GYY) 

90 CONTINUE 



CALCULATE THE FIRST-ORDER ALPHA AND BETA MATRICES 



ALPHA ( I , J)=( I./PI )*GNI J-DELTHE 
ALPHA ( J , I) = ( 1 ./PI ) *GN J I *D EL THE 
BETAtl, J)=(l./(2.*PI) )*GI J*DELTHE 
B ETA ( J » I ) = BET A ( I » J ) 

BET AX (I ,J) = ( l./(2.*PI) ) *GX I J^DEL THE 
B ETAX ( J , I )=( 1 ./( 2.*PI ) ) *3X J I *DEL THE 
BETAYt I , J) =( 1 ./< 2 .*PI ) )*GYI J*DELTHE 
BETAYt J , I ) = < 1 ./ ( 2.*PI ) )*GYJ I*DEL THE 
100 CONTINUE 

DO 120 1 = 1 , NPTS 

ALPHA ( I , I) =ALPHA( I , I KCMPLXt 1.0, 0.0) 

BETAX( I , I ) =BE T AX (1,1) +ANX t I ) *CMPLX ( 0 . 5 , 0 . 0 ) 
BETAYt I , I) =BETAY( I , I KANYt I )*CMPLX(0.5, 0.0) 
120 CONTINUE 



GENERATION OF THE FIRST-ORDER DISTRIBUTION FUNCTION, 
F ( I , 1 ) , BY INVERSION OF ALP HA ( I , J ) * F ( I , 1 ) = HH( 1,1) 



CALL COMAT ( 24, 1 , ALPHA, HH,DET, INDICA) 



THE HH (1,1) MATRIX IS REPLACED BY THE DISTRIBUTION 
FUNCTION, Ft 1,1) 

COMPUTATION OF THE FIRST-ORDER SCATTERING POTENTIAL 
FUNCTION, UISC(I), AND ITS X AND Y PARTIAL DERIVATIVES 
UISCX(I) AND UISCY(I), AND COMBINATION WITH THE 
INCIDENT POTENTIAL COUNTERPARTS TO COMPUTE THE TOTAL 
POTENTIAL FUNCTION AND ITS X AND Y DERIVATIVES 



DO 155 1=1, NPTS 
SUM 1= ( 0 . 0 , 0. 0) 

SUM2=( 0 .0,0.0) 

SUM3=( 0 .0,0.0) 

DO 150 J=1,NPTS 

SUM1 = SUM1 + F l J , 1 ) *BETA( I , J ) 

SUM2=SUM2+Fl J,1)*BETAX( I, J) 

SUM3=SUM3+F( J,1 ) ^BETAYt I, J) 

150 CONTINUE 

Ull I ) =S UM1 - ( CHY ( I ) / ( A*CH AH ) ) *C E X P ( CM PLX < 0 . 0 , A*X l I ) ) ) 
U 1 X ( I ) =SUM2-CMPLX( 0 . 0 , 1 . 0 ) *CHY ( I ) *CEXP t CMP LX l 0 . 0 , 
1A*X< I ) ) J/CHAH 

UlYt I ) =SUM3-SHY ( I ) *C EX P ( CM PLX ( 0 . 0 , A*X t I ) ) J/CHAH 
155 CONTINUE 



EVALUATION OF THE FREE SURFACE PRESSURE DISTRIBUTION 
SOURCE STRENGTH FUNCTION, FS(L) 



DELX=0.015625*PI/A 
WRITE (6,180) DE LX 
130 FORMAT ( 5X//5X , » DELX= ' , F 1 0. 5/// ) 
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SP=0.0 

CST=.5*DELTHE/PI 
I COUNT = 1 

DU 250 1 = 1 » NS PTS 
SUM1=( 0 .0,0.0) 

SUM2=( 0 .0,0.0) 

SUM3=(0. 0,0.0) 

SUM4=(0 .0,0.0) 

DO 245 1=1 , NPTS 
XV=SP 
Y V=0 . 0 
XC=X< I ) 

YC=Y( I ) 

ANX 1 = 0 . 0 
ANY 1=1.0 
ANX J=AN X ( I ) 

ANY J= AN Y ( I ) 

SHY 1= SHY ( 2 5 ) 

CHY I=CHY( 25) 

SHY J = SHY( I ) 

CHY J = CHY ( I ) 

IF (ABS(XV-XC) .LE.SMIN) GO TO 240 

CALL GREENS ( A , ANU , SH2 AH, SHAH , CHAH , COSAMU , S INAMU , AMU , 
1C OEFG , S HY I , CHY I , SHY J , CHY J , 2 5 , I , XV , YV , XC , YC , ANX I , ANY I , 
2ANXJ, ANY J,GI J ,GNI J,GNJI ,GXI J,GXJ I ,GYI J,GYJI , GYY ) 

GO TO 241 

240 CONTINUE 

CALL GREEN ( A , ANU , SH2 AH , SHA H , CH AH, AO , AA , BB , CC , DD , 2 5 , I , 
1XV,YV,XC,YC, ANX I , ANY I ,G I J , GN I J , GX I J , GY I J , GY Y ) 

241 CONTINUE 

SUM1 = SUMH-CST*GI J*F( 1,1) 

SUM2=SUM2*CST*SXIJ*F( I , 1) 

SUM3=SUM3 + CST*GY I J*F( 1,1) 

SUM4=SUM4+CST*GYY*F( 1,1) 

245 CONTINUE 

U1IS=(-1./A)*CEXP(CMPLXO.O, A*SP) ) 
U1ISX=CMPLX(0.0»-1.0)*CEXP{CMPLX(0.0»A*SP) ) 

U1 ISY = -TANH( A*H) *CEXP( CMPLX (0.0 , A*SP) ) 
U1ISYY=-A*CEXP(CMPLX(0.0,A*SP) ) 

U1SS=SUM1 
U 1 SSX = SUM2 
U 1SSY=SUM3 
U 1SSYY= SUM4 

F S ( L ) = ( 2 . *A/ { 3 . * ANU ) ) * ( U1 I S*U1 SS YY +U 1 I S Y Y*U 1 SS + U 1 S S* 
1U1SSYY-6.*U1 ISY*U1SSY-3.*U1SSY*U1SSY-4.*U1 I SX*UISSX 
2-2.*UlSSX*UlSSX) 

IF ( ICOUNT.EQ.O ) GO TO 248 
S P=F LO AT ( L + 1 ) *DE LX/ 2 . 

I C0UNT=0 
GO TO 250 
248 SP=-SP 

I COUNT = 1 
250 CONTINUE 



CALCULATION OF THE PARTICULAR SOLUTION PORTION OF THE 
SECOND-ORDER SCATTERING POTENTIAL AND ITS NORMAL 
DERIVATIVE, U2SC1U) AND U2SCNKI) 



CST=.5*DELX/PI 
DO 350 1=1, NPTS 
SUM1=( 0 .0, 0.0) 
SUM2=( 0 .0 ,0.0) 
SP=0. 0 
I COUNT = 1 

DO 345 L=1 , NSPTS 
XV=X( I ) 
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YV=Y( I ) 

XS = SP 
Y S = 0 . 0 
ANX 1= ANX ( I ) 

ANYI= AN Y ( I ) 

ANX J=0 . 0 
ANY J= 1 . 0 
SHY 1= SH Y ( I ) 

CHY I =CHY ( I ) 

SHY J=SHY( 25) 

CHYJ=CHY( 25) 

IF (ABS (XV-XS) .LE.SMIN) GO TO 340 

CALL GREENS ( A2 , ANU4 , SH2AH2 , SHA H2 , CHAH2 , COS AM2 , S I NAM2, 
1AMU4,C0EFG2,SHYI , CHY I , SHY J , CHY J , I , 25 , X V , YV , XS , Y S , 
2ANXI,ANYI, ANXJ,ANYJ,GI J,GNI J,GNJI,GXI J,GXJI ,GYI J.GYJI, 
3GYY ) 

GO TO 341 

340 CONTINUE 

CALL GREEN ( A2 , ANU4 , SH2 AH2 , SHAH2 , CHAH2 , A 02 , AA2 , BB2 , CC2 
1 , DD2, I , 25, XV, YV, XS ,YS, ANX I , ANY I , GI J , GNI J , GX I J ,GY I J , 
2GYY ) 

341 CONTINUE 
SUM1=SUM1+CST*GIJ£FS( L) 

SUM2=SUM2*CST*GNI J*FS< L) 

IF ( ICOUNT.EQ.O) GO TO 290 
SP=FL0AT(L+l)*DELX/2. 

I C0UNT=0 
GO TO 345 
290 SP=-SP 

I COUNT = 1 
345 CONTINUE 

U2SC1 ( I ) =SUM I 
U2SCNK I ) = SUM2 
350 CONTINUE 



CALCULATION OF THE HOMOGENEOUS SOLUTION PORTION 
OF THE SECOND-ORDER SCATTERING POTENTIAL, U 2SC0 ( I ) 



DO 500 1=1 , NPTS 
DG 500 J=1,I 
XV=X( I ) 

Y V = Y ( I ) 

XC=X( J ) 

YC=Y ( J ) 

ANX I = ANX ( I ) 

ANYI = ANY ( I ) 

ANXJ=ANX( J ) 

ANY J= ANY ( J ) 

I F ( ABS ( X ( I )-X( J) ) .LT.SMIN) GO TO 450 
SHY I=SHY ( I ) 

CHY I = CH Y ( I ) 

SHY J= SH Y ( J ) 

CHY J=CH Y ( J ) 

CALL GREENS ( A2 , ANU4 , SH2AH2 , SHAH2 , CHAH2 , COS AM2 , S I NAM2, 
1AMU4,C0EFG2, SHY I ,CHYi , SHY J , CHY J , I , J , XV , Y V , X C , YC , ANX I , 

2 ANY I ,ANXJ»ANYJ,GIJ ,GNIJ,GNJ I »GXIJ,GXJ I , GYI J ,GYJ I , GYY ) 

GO TO 490 
450 CONTINUE 

CALL GREEN I A2 , ANU4 , S H2 AH2 , SHAH2 , CHAH2 , A02 , AA 2 , BB 2 , CC2 
1 ,DD2 , I , J,XV, YV, XC , YC. , ANX I , ANY I , G I J , GN I J , GX I J, GY I J , 

2GYY J 

IF(I.EQ.J) GN J I =GN I J 
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I F { I.EQ.J ) GO TO 490 



CALL GREEN ( A2 , ANU4, S H2 AH2 , SHAH 2 , CHAH2 , A02 , AA 2 , BB 2 , CC2 
1,DD2, J, I,XC,YC,XV,YV,ANXJ,ANYJ,GI JEXT,GNJI , GXJI,GYJI , 
2GYY) 

490 CONTINUE 

ALPHA ( I , J)=(1./PI)*GYIJ*DELTHE 
ALPHA! J,I)=< 1 . / P I ) *GY J I *D EL THE 
BETA! I , J) =( 1./ ( 2 .*PI ) )*GI J*DELTHE 
BETA! J, I) =BETA( I , J ) 

500 CONTINUE 

DO 510 1 = 1, NPTS 

PK! I » 1 ) =2* ( PK ( I , 1 ) -U2SCN1 ( I )) 

510 CONTINUE 

DO 520 1 = 1 » N PT S 

ALPHA! I , I ) =ALPHA ( I , D+CMPLX! 1.0, 0.0) 

520 CONTINUE 

CALL COMAT ( 24 , 1 , AL PH A , PK , D E T , I ND I CA ) 

DO 540 I =1 , NPTS 
SUM=( 0.0, 0.0) 

DO 530 J = 1 , N PTS 
SUM=SUM+Fl!J,l ) £ BET A! I , J) 

530 CONTINUE 

U2SC0! I ) =SUM 
540 CONTINUE 



EVALUATION OF THE TOTAL SECOND-ORDER POTENTIAL, U2(I) 



DO 550 1=1 , NPTS 

U2( I )=U2SC 1! I ) +U2SC0! I )-(CHY( I ) / ( 2 . * A* ( SH AH**4) ) )* 
1CEXP(CMPLX(0.0,2.*A*X( I ) ) ) 

550 CONTINUE 

WRITE (6,560) 

560 FORMAT ( 5X// 9X , • I • , 1 5 X , • U 1 ( I ) ' , 24X , • U1 X ! I ) • , 2 5X , 
1'UIYII ) ' , 2 4X , ' U2 ( I ) •//) 

DO 580 1=1, NPTS 

WRITE (6,570) I , U 1 ( I ) , U IX ( I ) , U 1 Y ( I ) , U2 ( I ) 

570 FORMAT ( 5X , I 5 , 4 ( F 1 6 . 6 , F 14 . 6 ) ) 

580 CONTINUE 



EVALUATION OF THE FIRST-ORDER, SECOND-ORDER PERIODIC, 
AND SECOND-ORDER STEADY STATE FORCE COEFFICIENTS AND 
THE PERIODIC FORCE PHASE SHIFT ANGLES 



SUM1=(0. 0,0.0) 

SUM2=(0 .0,0.0) 

SUM3=( 0.0,0. 0) 

S UM4= ( 0 .0, 0.0) 

SUM5=0. 0 

SUM6=0 . 0 

DO 720 1=1, NPTS 

SUM1 = S UM1 +U 1 ( I)*ANX( I >*DELTHE 

SUM2=SUM2+U1 ( I ) * ANY ( I ) *DELTHE 

SUM3 = SUM3 + ( ( 6.*ANU*ANU*U2( I ) /A ) -UiX ( I ) *U 1 X ( I ) — U 1 Y ( I ) 
1 *U1 Y ! I ) J^ANX! I)*DELTHE 

SUM4=SUM4+ ( !6.*ANU*ANU*U2( I)/A)-UlX!i)*UlX(I)-UlY(I) 
1 *U1Y ( I ) )*ANY( I ) *D EL F HE 

SUM5 = SUM5H CABS (U1X! I)*U1X( I ) ) + CA3S (U1Y! I )*U1Y( I) ) 
l+ANU*ANU/A-*2-l .0) *ANX ( I ) *DELTHE 
SUM6=SUM6+ l CABS ( U1 X ( I)*U1X{ I) ) +CABS(U1Y( I)*U1Y(I ) ) 
l+ANU-ANU/A#-2-l .0) *AN Y ( I )*DELTHE 
720 CONTINUE 

C1(1)=CABS(SUM1*A) 

C1(2)=CABS(SUM2*A) 

C 2 ( 1 ) = CABS ( SUM3*A* A/ ( 4.-ANU ) ) 
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C2(2)=CABS(SUM4*A*A/(4.*ANU) ) 

C3(1)=SUM5*A*A/(4.*ANU) 

C3(2)=SUM6*A*A/( 4.* ANU) 

PHASE 1 ( 1 ) = AT AN2 ( A I MAG ( SUM 1 * A ) f REAL(SUM1*A) ) 

PHASE I (2)=ATAN2( AIMAGt SUM2*A) f REAL(SUM2*A) ) 

PHASE 2 ( 1 ) = AT AN2 ( A I MAG ( SUM3*A*A/ (4.*ANU) ) , RE AL ( SUM3*A*A 
1/ (4.*ANU) ) ) 

PHASE2 ( 2)=ATAN2( AIMAG(SUM4*A*A/ (4.*ANU) ) , RE AL ( SUM4* A*A 
1/ (4.*ANU) ) ) 

WRITE (6*730) C 1 ( 1 ) , C 1 ( 2) , C2 ( 1 ) , C2 ( 2) , C 3 ( 1 ) »C3( 2 ) , 

1 PHASE I ( 1) , PHASE1 { 2 ) * PH ASE2 ( 1 ) , PHASE2 ( 2) 

730 FORMAT ( 5X///5X , • Cl ( 1 ) = • , F8 . 5 , 5X , • Cl ( 2 ) = • , F 8 . 5/ / 5X , 

1 • C2( 1) =' » F8.5.5X, *C2(2)=' , F8.5//5X, *C3( 1 )=• ,F8.5, 5X, 

2* C3 ( 2 ) = ' *F8.5//5X» 'PHASEll 1) = ' , F 8. 5 * 5X , ' PHASEI(2) = ' » 
3F3.5//5X* ' PHASE 2 ( 1 ) = • ,F8.5 f 5X, ' PHASE2(2)=* ,F8.5//) 

STOP 

END 



THIS SUBROUTINE READS THE INPUT GEOMETRICAL DATA AND 
CALCULATES THE FIRST 200 ROOTS OF AMU*T AN ( AMU*H) + 

ANU = 0 AND AMU4*TAN( AMU4*H) + ANU4 = 0. IT ALSO 
GENERATES ARRAYS OF CERTAIN FUNCTIONS AND COEFFICIENTS 
USED IN GREEN AND GREENS AS WELL AS THE HH AND PK 
MATRICES 



SUBROUTINE GEODAT ( A , A2 * ANU , ANU4* SH2AH , SH2AH2 * SHAH* 
1SHAH2,CHAH,CHAH2,A0,AA,B8,CC,DD,A02,AA2 , B62 ,CC2, 0D2 ) 

COMPLEX HH ( 24, I ) , PK ( 24 * I ) 

DIMENSION X( 24) , Y{ 24) ,ANX( 24) , ANY( 24) ,CHY(25) ,CHY2<25) 
1,SHY(25) * SHY2(25) , COEFG ( 200 ) , C0EFG2 ( 200 ) , AMU (200) , 

2 AMU4( 200) » COS AM J ( 20C , 25) , C0SAM2 ( 200 , 25 ) , S I N AMU( 20 0,25) 
3, SI NA M2 (200*25) , SH2Y ( 24) , CH2Y ( 24 ) ,XX ( 24 ) 

CCMMON/CPX/HH, PK 
CCMMON/VAR/X* Y, ANX, ANY 

COMMON/ GSHY/CHY,CHY 2, SHY , SH Y 2 , S H2Y , CH2 Y 
CQMMQN/GMU/CQSAMU ,C0SAM2 , S I NAMU , S I NAM2 , AMU , AMU4 , C OE FG, 
IC0EFG2 

COMMON/ CON ST/ H * D,DELTHE,SMIN*NPTS*NSPTS 
COMMON/CPI/PI 

RE AD (5, 5) A,H,D, SMIN*NPTS*NSPTS 
5 FORMAT (4F10 . 5, 2 14) 

ANU=A*TANH(A*H) 

S H2AH=S I NH ( 2 . * A* H ) 

CHAH=COSH( A*H ) 

S H A H= S I N H ( A*H ) 

AO=T AN H( A*H) 

AA=1./CHAH**2 
E B=-H*SHAH/(CHAH**2) 

CC=-H*H*( 1 ,-SHAH**2)/ (3.*CHAH**3) 

DD=H** 3* SHAH*( 2 .*C HAH** 2+3. *( 1 . - SHAH** 2 ) ) / ( CHAH**4* 1 2 ) 

CELTHE=2.*PI/NPTS 

THETA=DELTHE/2 . 

DO 15 I=1*NPTS 
X ( I ) = C0 S ( T HET A ) 

Y ( I) = SIN( THETAJ-D 
ANX ( I ) =X( I ) 

ANY ( I ) = Y ( I ) +D 

SHY ( I ) = S I NH ( A* ( H + Y ( I ) ) ) 

CHY( I )=COSh( A*( H+Y ( I ) ) ) 

S H2Y ( I)=SINH(2.*A*(H+Y( I ) ) ) 

CH2Y( I )=C0SH(2.*A*(H+Y( I ) ) ) „ , 

HH ( I , 1 ) = ( 2 ,/CHAH) *CMPLX( ANY (I)*SHY(I),ANX(I)*CHY(I))* 
1CEXP( CMPLX(0.0 , A*X l I ) ) ) 

PK( 1 , 1 ) =( ANY ( I ) * S H 2 Y ( I ) +CMP LX ( 0 . 0 , 1 . 0 ) * ANX ( I ) *CH2 Y ( I ) ) 
1*CEXP(CMPLX(0.0, 2 .0*A*X(I )) )/SHAH**4 
THETA=THETA+DELTHE 
15 CONTINUE 

SHY ( 25 ) =SH AH 
CHY(25)=ChAH 
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B=ANU*H 
DO 50 K=l, 200 
XX( 1J=PI^K 
DO 25 1=1,20 

I 1 = 1 

Y Y=XX ( I ) 

XX ( I+1)=XX( 1 ) - ATAN2(B ,YY) 

I F ( AB S ( (XX(I+1)-XX(I))/XX(I+1)) .LT..0001) GO TO 30 
25 CONTINUE 

IF ( II.GE.20) GO TO 27 
GO TO 30 

27 WRITE (6,28) 

28 FORMAT (5X,'AMU ROOT DOES NOT CONVERGE'//) 

30 CONTINUE 

AMU(K) =XX ( I I ) / H 

C0EFG(K)=2.*PI*( AMU(K) *AMU( K) +ANU*ANU)/( ANU*AMU( K) -H* 
1AMU(K)*(AMU(K)*AMU(K)+ANU*ANU) ) 

N PTS1 = N PTS +- 1 
DO 40 I = 1 , NP TS 1 
IF (I.EC.NPTS1) GO TO 35 
COSAMU ( K , I )=COS(AMU(K)*(H+Y ( I ) ) ) 

SINAMU< K, I) = SIN( AMU(K )*(H + Y( I ) ) ) 

GO TO 4C 

35 COSAMU ( K. , I )=COS ( AMU ( K ) *H) 

S I N AMU ( K , I )=SiN(AMU(K)*H) 

40 CONTINUE 
50 CONTINUE 

ANU4=4.*ANU 
B2=ANU4*H 
DO 70 K= 1 , 200 
XX( 1 )=PI*K 
DO 55 1=1,20 

I I = I 

Y Y=XX ( I ) 

XX ( 1+1) =XX< 1 ) -ATAN2(B2 ,YY ) 

I F ( AB S ( (XX( I + 1 ) - XX ( I ) ) /XX ( I +1 i ) .LT..0001) GO TO 60 
55 CONTINUE 

IF ( II.GE.20) GO TO 57 
GO TO 6 C 

57 WRITE (6,58) 

58 FORMAT (5X,'AMU4 ROOT DOES NOT CONVERGE'//) 

60 CONTINUE 

AMU4(K) =XX( I I )/H 

CGEFG2(K)=2.*PI*(AMU4(K) ! * : AMU4(K) +ANU4* ANU4 ) /( ANU4* 

1 AMU4 ( K ) -H*AMU4( K )*( AMU4 ( K ) * AMU4 ( K) +ANU4* ANU 4) ) 
NPTS1=NPTS+1 
DO 65 I = 1 , NP TS 1 
IF (I.EQ.NPTS1) GO TO 68 
C0SAM2 ( K , I )=COS( AMU4 ( K ) * ( H+ Y ( I ) ) ) 

SINAM2 (K, I ) = S I N ( AMU4 (K)*(H + Y( I ) ) ) 

GO TO 65 

68 C0SAM2(K,I ) =COS ( AMU4( K ) *H) 

S I NAM2 ( K, I)=SIN(AMU4(K)*H) 

65 CONTINUE 
70 CONTINUE 
XX ( 1 ) =ANU4 
DO 80 1=1,20 
I I = I 

Y 2=XX ( I ) 

XX ( 1+1 ) =4. *ANU/TANH(Y2*H) 

IF (A3S(lXX( I+1)-XX( I))/XX( 1 + 1) ) .LT. 0.3001) GO TO 85 
80 CONTINUE 

WRITE (6,32) 

82 FORMAT (5X,'A2 ROOT DOES NOT CONVERGE'//) 

85 CONTINUE 
A2=XX ( I I ) 

SH2AH2=SINH( 2.*A2*H) 

CHAH2=C0SH( A2*H) 

SHAH2=S INH( A2*H) 

A02=TANH( A2*H) 

AA2=1./CHAH2**2 
BB2=-H*SHAH2/ ( CHAH2**2 ) 
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CC2=-H*H*{ l.-SHAH2**2) /(3.*CHAH2**3) 
DD2=(H**3)*SHAH2*( 2. *CHAH2**2+3 . *( l.-SHAH2**2))/ 

1 ( (CHAH2**4)*12. ) 

DO 90 I = 1 , NPTS 

SHY2( I)=SINH(A2*(H+Y( I)) ) 

C HY 2 ( I )=C0SH(A2*(H+Y( I) ) ) 

90 CONTINUE 

WRITE (6,95) A,A2 ,D , H , ANU , A NU4 , NPTS , SMI N , NS PTS 
95 FORMAT ( 5X///3 X , ' A = • , F 1 0 . 5 , 4X , ' A 2 = ' , F 10 . 5 , 4X , ' D= ' , 
1F10.5,4X,'H=',F10.5,4X, ' ANU=' , F 1 0 . 5 , 4X , ’ ANU4= • ,FI0.5/ 
24X , 1 NPT S= ' , I3,4X , ' SMIN=' , FI 0.5 ,4X, , NSPTS = * , 13/// ) 
RETURN 
END 



THIS SUBROUTINE CALCULATES G AND ITS DERIVATIVES, 

GX , GY , GN , AND GYY, 3Y USE OF THE INTEGRAL FORM FOR THE 
CASE (X(I) - X(J)) LESS THAN SMIN 



SUBROUTINE GREEN ( A , ANU , S H2 AH, SH AH ,CHAH , AO , AA , SB , CC , DD 
1, I , J,X, Y,XI ,ETA, ANX,ANY,G,GN,GX, GY, GYY) 



COMPLEX G, GN 
COMPLEX GX ,GY » GY Y 

DIMENSION TEST! 200) , TESTT( 100) , TE STTT ( 100) ,SUMOX( 15) , 

1 SUMOY (15) ,NNN(15) ,TESTYY(200) 

C CMMON/ CON ST/H,D,DELTHE, SMIN, NPTS, NS PTS 
COMMON/CPI/PI 

PI (Y »ETA, XX, AMU) = C OSH (AMU* ( Y + H ) )*COSH(AMU*( ETA + H) ) * 
1C0S( AMU* XX ) /(COSH ( AMU *H)**2 ) 

P2X( Y,ETA,X,XI , AMU)=-AMU*COSH( A MU* (ETA+H ) ) *SIN( AMU* 

1 (X-X I ) )*COSH( AMU* (Y+H) ) / ( COSH ( AMU*H ) **2 ) 

P2Y(Y,ETA,X,XI ,AMU)=AMU*COSH( AMU* (ETA+H ) ) *S INH( AMU* 

1 (Y + H) ) *COS (AMU* (X-X I ) ) / ( C OS H( A MU*H ) **2 ) 

P 2YY (Y,ETA,X*XI , AMU ) = A MU*AMU*C OS H ( AM J* ( Y +H ) )*CCSH( AMU* 
1 ( ETA+H) ) *COS( AMU*( X-X I ) ) / ( COSH ( AMU*H ) **2 ) 

Q1 ( Y ,ETA, XX, AMJ) =-Pl ( Y, ETA, XX , AMU) *( AMU-A)/ ( AMU* 

IT ANH( AMU*H)-ANU) 

Q2X(Y,ETA,X,XI,AMU)=-P2X(Y,ETA,X,XI, AMU) *( AMU-A) /(AMU* 
1 T ANH (AMU*H)-ANU) 

Q2Y(Y , ETA, X, XI, AMU ) =-P2Y ( Y , ETA , X , X I , AMU ) * ( A MU- A) /(AMU* 
IT ANH( AMU*H)-ANU) 

Q2YY( Y, ETA,X,XI , AMU)=-P2YY (Y,ETA,X,XI , AMU ) * ( AMU- A ) / 

1 ( AMU*T ANH( AMU*H ) -ANU) 

Q 10 ( Y ,ETA,XX)=-COSH(A*(Y+H) ) *COS H ( A* ( ET A + H) )*COS(A*XX) 
1*A/(C0SH(A*H)**2*( ( A* A -ANU* ANU) *H+ANU) ) 

Q2CX(Y , ET A , X , X I ) =A*COSH ( A* ( ETA+H ) ) *S I N ( A* ( X -X I ) )*COSH 
1 ( A*( Y+H ) ) *A/(COSH( A*H)**2*( ( A* A- ANU* AN J ) *H+ ANU) ) 

Q 20Y (Y,ETA»X»XI )=-A*COSH(A*( ETA + H) ) * S I NH ( A* ( Y +H ) ) *C OS 
1 ( A* ( X-X I ) )*A/(C0SH(A*H)**2*( ( A* A- ANU* ANU ) *H + ANU) ) 

Q20YY (Y , E T A , X , X I ) =-A*A*COSH ( A* ( Y + H) )*C3SH(A*( ETA + H))* 
1C0S(A*( X-XI ) )*A/(C0SH(A*H)**2*( ( A* A- ANU* ANU ) *H + ANU ) ) 

Q IS (Y, ETA, XX, A MU) =-Pl ( Y , ETA , XX , AMU) / ( A0+ AMU*H* ( AA +BB* 

1 ( AMU-A) +CC*( AMU-A)**2+DD*( AMU-A )**3) ) 

Q2XS ( Y , ETA,X,XI,AMU)-=-P2X(Y,ETA,X,Xl,AMU)/(A0+AMU*H* 

1 (AA+BB* (AMU-A) +CC*( AMU- A ) **2 +DD* ( AMU- A ) * * 3 ) ) 

Q2YS( Y, ETA ,X, XI , AMU ) =-P2Y ( Y ,ETA , X , XI , AMU ) / ( A0+AMU*H* 

1 (AA+BB* (AMU-A) +CC*( AMU-A) ** 2 +DD* ( AMU-A) **3) ) 

Q2YYS( Y , ETA, X.XI , AMU) =-P2YY (Y,ETA,X,XI, AMU ) / ( A0+ AMU*H* 
1 (AA+BB* (AMU-A )+;C*( AMU-A )** 2 +DD* ( AMU-A ) **3 ) ) 

FUNK Y, ETA, XX, AMU) =EX P ( -A MU*H ) * S I NH( AMU* ET A ) * S I NH 
1 ( AMU*Y) *COS( AMU* XX) / ( AMU*COSH( AMU*H ) ) 

FUN2( Y, ETA, XX) = 4 . * 3 . 1 4 1 59*C OSH ( A* ( Y+H ) ) *CCSH( A* 

1( ETA+H) )*C0S(A*XX)/(2.*A*H+SH2AH) 

FUNR(X, Y,XI ,ETA, aNX,ANY)=( (X-X I ) *ANX+ ( Y +ET A ) *ANY ) / 
l ( (X-XI )**2 + (Y + ETA)**2) 

FUN3X ( Y , E T A, X , X I , AMU ) =EXP ( - AMU*H ) *S I NH ( AMU* ETA) *S INH 
1 ( AMU*Y ) * S I N ( AMU* (X-X I) ) /C OS H ( AMU *H ) 

FUN3Y ( Y , ETA, X, XI ,AMU) =-EXP ( -AMU*H ) *S I NH ( AMU*ETA ) *COSH 
1 (AMU*Y)*COS( AMU* (X-XI) )/COSH(AMU*H) 

FUN3YY( Y,ETA,X,XI , AMU ) =-E X P ( - AMU*H) *S I NH ( AMU* ET A ) * 
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IS INH( AMU*Y ) *AMU*AMU*COS ( AMU* (X -XI ) )/COSH(AMU*H) 
FUN4( Y, ETA,XX,ANX,ANY , SIGN) =4. *3. 14159* A*COSH ( A* 

L ( ETA+H) ) *( ANY*SINH( A*< Y+H I ) *COS( A*XX) -ANX*COSH( A* 
2 ( Y+H ) ) *SIN(A*XX)*SIGN)/(2.*A*H+SH2AH) 



EVALUATION OF THE FINITE INTEGRAL IN G TO DETERMINE 
THE SIZE OF THE SUBDIVISIONS 



XX=ABS ( X-X I ) 

IF(X.LT.XI) SIS N =-1.0 
IF(X.GT.XI) S I GN = 1 . 0 
IF(X.EQ.XI) S I GN = 0 . 0 
DO 50 N=l,15 
DELMU=2*A/(o*N+3 ) 

AMU=0 . 0 
SUM=0 . 0 

F0=(Q1( Y, ETA, XX, AMU)-Q10(Y , ETA, XX) )/( AMU- A) 

LL=6*N+3 
DO 40 NN= 1 » LL 

IF( ABS( AMU-A) ,LT . .00001 ) GO TO 10 
IF(ABS( AMU+DELMU-A) .LT. .00001) GO TO 10 
IF(A6S( AMU+DELMU/3. -A) .LT . .00001 ) GO TO 10 
IF(ASS( AMU+2.*DELMU/3.-A) .LT. .00001) GO TO 10 
F 1= ( Q 1 ( Y, ETA, XX, AMU+DE LMU/3. ) -Q 1 0 ( Y ,E TA , XX ) )/< AMU+ 
1DELMU/3.-A) 

F 2= ( Q1 ( Y,ETA,XX,AMU+2.*DELMU/3. ) -Q10 ( Y , ETA, XX ) ) / 

1( AMU+2.*DELMU/3.-A) 

F3=(Q1I Y, ETA, XX, AMU+DE LMU J -Q10C Y , ETA ,XX ) )/( AMU+DE LMU 
1 — A J 

GO TO 30 

10 F0=(Q1(Y,ETA,XX, AMU+DE LMU ) -Q10 ( Y , ETA , XX ) ) /( AMU+DE LMU 
1- A ) 

GO TO 40 

30 SUM= (DELMU/8. ) * ( FO +3 .*F 1 +3. *F2 +F3 ) +SU M 
F 0=F 3 

40 AMU= AMU+DE LMU 
TEST ( N ) =SUM 
IF(N-l) 50,50,45 
45 MN=N-1 

MM=6*MN+3 

IF(ABS( (TEST(N)-TEST(N-l) )/TEST(N) ) .LT..010) GO TO 60 
50 CONTINUE 
60 CONTINUE 
PV1F=2.*SUM 



EVALUATION OF THE FINITE INTEGRAL IN GN USING 
2* A/MM SUBDIVISION SIZE 



DELMU=2*A/MM 
AMU=0. 0 
SUMX=0 . 0 
SUMY=0 . 0 

F0=(Q2X(Y,ETA,X,XI ,AMU)-Q20X(Y ,ETA,X,XI ) )/( AMU-A) 
Y0=(Q2Y(Y,ETA,X,XI , AMU ) -3 20Y ( Y , ET A ,X , X I ) )/( AMU-A) 

DO 80 NN= l , MM 

IF(ABS( AMU-A) .LT. .00001) GO TO 70 

IF( ABS( AMU+DELMU-A) .LT. .00001) GO TO 70 

IF(A6S< AMU+DELMU/3 .-A ) .LT. . 00001 ) GO TO 70 

I F ( ABS( AMU+2.*DELMU/3.-A) .LT. .00001) GO TO 70 

F1=(Q2X(Y,ETA,X,X I , AMU +DELMU/ 3 . ) -Q20X ( Y ,ETA,X,XI))/ 

1 ( AMU+DELMU/3. -A) 

F2=(Q2X( Y ,ETA,X,XI , AMU+2. *0E LMU/3 . ) -Q20X ( Y , ETA,X, XI ) )/ 
1 ( AMU+2. *DEL MU/3. -A) 

F3=(Q2X(Y,ETA,X,XI , AMU + DELMU ) -Q2 OX (Y , ET A , X , X I ) )/ 

1( AMU+DELMU-A) 

Y1=(Q2Y (Y,ETA,X,XI , AMU+DE LMU/3 . ) -Q20Y ( Y ,ETA,X,XI ) )/ 

1 ( AMU+DELMU/3. -A) 

Y2=(Q2Y (Y , ETA,X»X I »AMU+2.*DELMU/3 . ) -Q20Y ( Y , ETA,X,XI ) )/ 
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1 ( AMU+2 . #DELMU/3.-A) 

Y3=(Q2Y ( Y, ETA,X,X [ , AMU+DELMU ) -rQ20Y (Y,ETA,X,XI))/ 
1 ( AMU+DELMU-A ) 

GO TO 75 
70 CONTINUE 

E 0= ( Q2X ( Y , ET A , X , X I , AMU +DE L MU ) - Q2 CX ( Y , ET A , X , X I ) )/ 
1 ( AMU+DELMU-A ) 

Y0= ( Q2Y (Y, ETA,X,XI, AMU+DELMU )-Q2CY(Y,ET A, X,X I ) )/ 
1 ( AMU+DELMU-A) 

GO TO 30 
75 CONTINUE 

SUMX= SUMX + (DELMU/8. ) * ( FO +3 . * F 1 +3 . *F2 + F 3 ) 

SUMY=SUMY+ (DELMU/8 . )* { YO +3 . *Y 1 +3 . *Y2 +Y3 ) 

F0 = F3 
YO = Y 3 

80 AMU=AMU +DELMU 
PV2FX=2 .*SUMX 
P V2FY=2 . *SUMY 



EVALUATION OF THE FINITE INTEGRAL IN GY Y USING 
2* A/MM SUBDIVISION SIZE 



DELMU=2*A/MM 
AMU=0 . 0 
SUMYY = 0.0 

YY0=(Q2YY( Y,ETA,X,XI , AMU) -Q 2 OY Y ( Y , ET A , X , X I ) J/IAMU-A) 

DO 350 NN= 1 , MM 

I F { ABS ( AMU- A ) . LT .0.00001) GO TO 320 
IF (ABS( AMU + DELMU-A) .LT. 0.00001 ) GO TO 3 20 
IF ( ABS ( AMU + DELMU/3. -A) .LT. 0.00001 ) GO TO 320 
IF (ABS(AMU+2.*DELMU/3.-A) .LT. 0.00001) GO TO 320 
YY1=(Q2YY(Y, ETA, X, XI, AMU+DELMU/3. )-Q2QYY(Y, ETA, X, XI ) )/ 
1 (AMU+DELMU/3. -A) 

YY2=(32YY{ Y, ETA,X ,XI , AMU+2 . *DE LMU/3. ) *3 2 OYY ( Y , E T A , X , 

IX I ) ) /( A MU + 2. +DELMU/ 3. -A) 

YY3=(Q2YY(Y,ETA,X,XI,AMU+DELMU)-G20YY(Y ,ETA,X,XI ) )/ 

1 (AMU+DELMU-A) 

GO TO 325 
320 CONTINUE 

YY0=(Q2YY(Y, ETA,X,XI , AMU + DELMU ) -Q20YY ( Y ,ETA,X» XI ) )/ 

1 ( AMU+DELMU-A) 

GO TO 350 
325 CONTINUE 

SUMYY =SUMYY+ (DELMU/8. ) * ( Y YO+3 . * YY 1 +3 . * Y Y 2+Y Y3 ) 

Y Y0=YY3 

350 AMU=AMU+DELMU 
PV2FYY=2.* SUMYY 



EVALUATION OF THE INFINITE INTEGRAL IN G, GN , AND GYY 
SIMULTANEOUSLY 



AMU=2* A 
DELMUO=DELMU 

F0=Q1 ( Y , ETA, XX, AMU) /( AMU- A) 

F0X=32X (Y,ETA, X, XI , AMU ) / ( AMU- A ) 

F0Y=Q2Y (Y »ETA»X,XI , AMU) /( AMU-A) 

F0YY=Q2YY( Y,ETA,X,XI ,AMU) /( AMU-A) 

SUM=0 . 0 
SUMX=0 . 0 
SUMY=0 . 0 
SUMYY=0 .0 
DO 100 NN= 1,200 
DO 95 N = 1 , 20 

F 1=Q1 ( Y , ETA, XX, AMU+DELMU/3. )/( AMU +DELMU / 3 .- A ) 

F 1X=Q2X ( Y , ETA ,X , XI , AMU+DELMU/3 . ) / (AMU+DELMU/ A .-A) 
F1Y=Q2Y(Y»ETA»X,XI , AMU+DELMU/3 . ) /(AMU+DELMU/3. -A) 
F1YY=Q2YY( Y,ETA ,X ,XI , AMU+DELMU/3. ) / ( AMJ + DE L MU/3 . - A ) 
F2=Q1 ( Y ,£TA, XX, AMU+2. *DELMU/3. ) / ( AMU+2 . * DEL MU/ 3 . -A ) 
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F2X=Q2X(Y,ETA,X,XI,AMU+2.*0ELMU/3. ) / ( AMU + 2 . *DELMU/ 3 . 
1-AJ 

F2Y=Q2Y(Y , ETA, X, XI , AMU + 2 . *DELMU/3 . ) / ( AMU + 2 . *DELMU/ 3 . 

I — A) 

F2YY=Q2YY(Y, ETA, X,XI ,AMLH-2.*DELMU/3 . ) / ( AMU + 2. *DELMU/ 3. 
1 - A) 

F3=Q1(Y»ETA,XX,AMU+ OELMU ) / { AMU+DELMJ-A J 
F3X=Q2X < Y , ETA , X , X I , AMU+DE LMU) / ( AMU+DELMU- A ) 
F3Y=Q2Y(Y,ETA,X,XI,AMU+DELMU)/(AMU+DELMU-A) 
F3YY=Q2YY(Y, ETA, X, XI, AMU+DE LMU ) / ( AMU+DE LMU- A) 

SUM= (OELMU/3. ) * ( FO +3 . *F 1 +3. *F 2 +F3 ) +SUM 
SUMX= ( DEL MU/ 8 . ) * ( FOX + 3 . * FIX + 3 . *F 2X + F 3 X ) + SUMX 

SUMY= (DELMU/3. ) * ( F0Y+3 . *FlY+3 . *F2Y+F3 Y ) +SUMY 

SUMYY=SUMYY+ (OELMU/8. ) *(FOYY+3 . * F 1 YY + 3 . * F2Y Y + F3 YY ) 

FO = F3 
F0X=F3X 
FO Y=F3Y 
F 0YY=F3 YY 

ASM=EXP(AMU*(ETA+Y) )/AMU 
IFIASM.LT. . OOOi ) GO TO 86 

95 AMU= AMU + DE LMU 
86 TEST! NN ) = SUM 

TESTT(NN) =SUMX 
TESTTT { NN) =SUMY 
TESTYY ( NN ) =SUMYY 
IF(NN-l) 100,100,97 

97 CONTINUE 

IFIASM.LT. .00010) GO TO 105 
I F ( ABS ( SUM) .LT. .0001) GO TO 98 

IF ( A 3 S ( (TEST (NN)-TEST (NN-1 ) )/TEST (NN) ) .GT..001) 

1G0 TO 100 

98 I F ( ABS C SUMX) .LT. .0001) GO TO 99 

IF ( ABS ( (TESTT(NN)-TcSTT(NN-l) ) / TESTT ( NN ) ) . GT. . 00 1 ) 

1 GO TO 100 

99 IF(A6S( SUMY) .LT. .0001) GO TO 93 

I F ( ABS ( ( TESTTT ( NN) -TESTTT ( NN-1) )/ TESTTT ( NN )). GT. . 001 ) 
1G0 TO 100 

93 IF(ABS( SUMYY) .LT. 0.0001) GO TO 96 

IF ( ABS (( TESTYY (NN) -TESTY Y( NN-1 ) ) / TE ST Y Y ( NN ) ) .GT.. 001) 
1 GO TO 100 

96 CONTINUE 
GO TO 105 

100 CONTINUE 

102 WR I TE(6,103) 

103 FORMAT (3X3 4H INFINITE INTEGRAL DID NGT CONVERGE) 

105 CONTINUE 

PV1I=2.*SUM 
PV2 IX = 2 .*SUMX 
PV2IY=2 .*SUMY 
PV2IYY=2.*SUMYY 
DELMU=. 2/H 

IF (J.EQ.25) GO TO 240 

AMU=0 . 0 

SUM=0 . 0 

SUMX=0 . 0 

SUMY=0 . 0 

SUMYY=0 .0 

F0=0 .0 



F0X=FUN3X(Y,ETA,X,XI, AMU ) 
F0Y=FUN3Y(Y,ETA,X,XI ,AMU) 
F0YY=FUN3YY(Y,ETA,X,XI ,AMU) 

DO 200 NN= 1 ,100 
I F ( AMU*H. GT . 5 . ) DELMU=.l 
DO 195 N=1 ,20 

F 1X = FUN3X( Y ,ETa, X,XI , AMU + DE LMU/ 3 . ) 
F2X=FUN3X( Y,ETA, X, XI , AMU+2. ^DELMU/3. ) 
F3X = FUN3X( Y ,ETA,X,XI , AMU + DELMU ) 
F1Y=FUN3Y( Y, ETA, X, XI , AMU+DE LMU/3. ) 

F 2Y= FUN 3 Y ( Y , ET A , X , X I , AMU+ 2 . *DE LMU/3. ) 
F3Y=FUN3Y( Y,ETA,X,XI , AMU+DELMU) 

F l YY = FUN3YY { Y , ET A , X , X I , AMU + DEL MU/ 3 . ) 
F2YY=FUN3YY (Y,ET A, X,X I , AMU+2. *DE LMU/3. ) 
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F3YY=FUN3YY( Y,ETA,X,XI , AMU+DELMU) 

I F ( AMU + DE LMU .LT . .001) GO TO 120 
F 1 = FUN 1 (Y, ETA, XX, AMU+DE LMU/ 3 . ) 

F2=FUN1 (Y,ETA,XX,AMU+2.*DELMU/3. 1 
F3=FUN1 (Y , ETA, XX , AMU+DELMU ) 

GO TO 130 
120 CONTINUE 

Fl = EXPl-( AMU + DELMU/3 . )*H)*CQS( ( AMU+DELMU / 3 . )*XX)* 

1( AMU+DELMU/3. ) *Y*E T A/ COSH ( ( AMU+DELMU/3. )*H) 
F2=EXP(-(AMU+2.*DELMU/3.}*H)*CGS( ( AMU+2 . *DE LMU/ 3 . }*XX) 
l*(AMU+2.*DELMU/3. ) *Y*ETA/COSH( ( AMU + 2 .*DE LMU/3 . ) *H ) 
F3=EXP(-( AMU+DE LMU) *H)*COS( (AMU+DELMU)-XX)* ( AMU+DE LMU } 
1 *Y*ETA / CO SH ( ( AMU+DE LMU ) *H ) 

130 CONTINUE 

SUM= ( DELMU/8. )* ( FO +3.*F1 +3 .*F2 +F3 } + SUM 

SUMX= ( DELMU/8 . )*( FOX +3 . * F IX +3 .-F2X + F 3 X } + SUMX 

SUMY= (DELMU/8. )*(F0Y+3.*F1 Y +3. *F 2Y+F 3Y ) + SUMY 

SUMYY=SUMYY+( DELMU/8 .)*(F0YY+3.*FIYY+3.*F2YY+F3YY) 

F0=F3 

F0X=F3X 

F0Y=F3Y 

F0YY=F3 YY 

ASM=EX P( AMU* ( ETA + Y ) ) 

IFtASM.LT. .0001 ) GO TO 196 



195 

196 



199 



A MU= AMU+DE LMU 
TEST(NN)=SUM 
TESTT(NN) =SUMX 
TESTTT ( NN ) =SUMY 
TESTYY ( NN) =SUMYY 
IF(NN-l) 200,200, 
CONT INUE 

IF(ASM.LT. .00010) 



199 

GO 



GT. .0010) 



TO 205 

IF(ABS( SUM) .LT. .0001 ) GO TO 2C6 
IFUBSl ( TEST(NN)-TEST INN-1) ) / TEST ( NN ) ) 

1 GO TO 200 

206 I F ( ABS ( SUMX) .LT. .0001) GO TO 207 

IF ( ABS ( ( TESTT(NN)-TESTT( NN- 1 ) ) / TE STT ( NN ) ) . G T . . 00 1 0 ) 
1G0 TO 200 

207 IF (ABS< SUMY) .LT. .0001 ) GO TO 204 
IF (ABS ( ( TESTTT ( NN)- TESTTT (NN-1 ) ) /TESTTT (NN) ) .GT 

1 GO TO 200 

204 IF( ABS( SUMYY) .LT. 0.0001) GO TO 203 

IF (ABS (( TESTYY ( NN )-TESTYY( NN-1 ) ) / ( TEST Y Y ( NN) ) ) 

1GT. .001 ) GO TO 200 



0010 ) 



208 

200 

202 

205 



CONVERGENCE) 



CONTINUE 
GO TO 205 
CONT INUE 
WRITE(6 ,202) 

FORMAT ( 3X14HN0 
CONT INUE 
GINF=-2*SUM 
GX I NF=2 . *SUMX 
GYINF=2.*SUMY 
GYYINF=2.*SUMYY 
GNSING= .5 

IF (I.EQ.25) GO TO 218 
I F( I.EQ . J) GO TO 220 
A I J J- I - J 

THETA1= ABS( AIJJ) *DELTHE 

AINJJ=I +NPTS- J 

THETA2 = ASSl A INJJ UDELTHE 

A JN 1 1 = I - J- NPT S 

THETA3 = A3S( AJNI I ) *DEL THE 

THETA=THETA1 

IF (THETA2.LT. THETA) THETA=THETA2 
IF(THETA3.LT. THETA) THET A = THET A 3 
I F ( THE TA.GT..15) GO TO 218 

GSING=DELTHt*ALOG( 2. ) + 2 . * ( -DEL THE/ 2. + ( D EL THE/ 4. + 
1THETA/2 . ) * ALOG ( THETA/2. +DELTHE/4. ) - ( THE T A/ 2 .-DEL THE / 4 . 
2 ) *ALOG ( THETA/2. -DEL THE/ 4. ) - ( DE LTHE/4 . + TH ET A/ 2 . ) ** 3/ 1 3 . 
3+ (THET A/ 2. -DEL THE /4. )**3/18.-(THETA/2.+DELTHE/4.)**5/ 
4( 180 .*5 . ) + ( THET A/ 2. -DEL THE/ 4. ) **5/ ( 180 . *5. ) -( THET A/2 . + 
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5DELTHE/4. ) ** 7/ ( 2 83 5. *7. ) + ( THETA/ 2. -DELTHE/4. ) -*7/ 
6(2835. *7.) ) 

GXSING= ( X-XI )/( ( X-XI }**2-HY-ETA)**2) 

GYSING=(Y-ETA)/( (X-XI ) **2+( Y-ETA ) **2) 

GO TO 230 

218 GSING=DELTHE*ALOG(SQRT( (X-XI } **2 + ( Y-ET A ) **2 ) ) 

GXSING= ( X-XI )/( (X-XI )**2*( Y-ETA )**2) 

G YS I NG= (Y-ETA) /( (X-XI ) **2 + ( V-ET A ) **2 ) 

GO TO 230 
220 CONTINUE 

G SI NG=D E LTHE* ALOG (DELTHE/2. ) -D ELTH E- ( DE L THE ** 3 ) / 

1( l8.*16.)-(DELTHE**5)/( 180 . *5. * 2 56 . ) - ( D E LTH E**7 ) / 

2( 2835. -7. *256. *16. ) 

GXSING=GNS ING*ANX 
GYSING=GNSING*ANY 
230 CONTINUE 

G =GSI NG/DcLTHE-ALOG(S QRT ( (X-XI )**2 + (Y + ET A )**2 ) ) + PV1F 
H-PV1I+GINF+CMPLX (0. , -1. )*FUN2(Y , ETA, XX) 

G X=GX S I NG- FUNP ( X , Y , X I , E T A , 1 . 0 , 0 . 0 ) +P V 2F X * PV 2 1 X <-GX I NF + 
1FUN4( Y , ETA ,XX, I . 0 , 0 . 0 , S IGN ) *CMPLX (0 . 0 , - 1 . 0 ) 
GY=GYSING-FUNR( X ,Y ,XI , ETA , 0 .0 , 1 . 0 ) +PV2FY +PV 21 Y+GY I NF+ 
1FUN4(Y , ETA, XX, 0. 0,1.0, SIGN) *CMPLX (0.0, -1.0) 

IF (I.LE.NPTS) 30 TO 235 

FURYY= 1 ./SORT ( ( ( X-XI ) **2+ ( Y-ET A ) **2 ) **3 ) - ( ( Y-ETA ) **2 + 
1( Y-ETA) ) /SORT ( ( (X-XI ) **2 H Y-ET A ) ** 2 ) **5 ) 

FUNYY= 1 . / SQRT ( ( (X-XI ) **2+( Y+ETA ) **2) **3 ) -( ( Y+ETA)**2+ 

1 ( Y + ETA) )/SQRT( l (X-XI )**2MY+ETA)**2)**5) 
GYY=FURYY*FUNYY+PV2FYY+PV2I YY+GYY INF+CMPLX( 0.0, -I . 0)* 
1FUN2(Y , ETA , XX } *A*A 
235 CONTINUE 

IF (I.EQ.25) GO TO 250 

GN=GNS I NG-FUNR ( X , Y , X I , ETA , A NX , ANY ) + ( PV2F X + P V2 IX+GX I NF ) 
l*ANX+{ PV2FY+PV2IY+GYINF)*ANY+FUN4(Y,ETA,XX, ANX, ANY, 
2SIGN)*CMPLX(0.0,-1.0) 

GO TO 250 
240 CONTINUE 

G=- (PV1F+PV1I )+CMPLX(0. ,-l. )*FUN2(Y,ETA,XX) 

GN=-( PV2FX+PV2IX ) * ANX- ( PV2F Y+PV2 I Y ) *ANY + FUN4( Y , ET A , XX, 
1 ANX , A NY , S I GN ) *C M PL X ( 0 . 0 , - 1 . 0 ) 

250 CONTINUE 
RETURN 
END 



THIS SUBROUTINE CALCULATES G AND ITS DERIVATIVES, 

GX , GY , GN , AND 3YY , BY USE OF THE SERIES FORM FOR THE 
CASE ( X ( I ) - X ( J ) ) GREATER THAN SMIN 



SUBROUTINE GREENS ( A , ANU, SH2 AH , SHAH ,CHAH , COSA MU , S I N AMU 
1 , AMU,COEFG,SHYI ,CHYI, SHYJ,CHYJ, I , J , XV , Y V , X C , YC , AN X I , 

2 ANY I , ANX J , ANY J , G I J , GN I J , GNJ I , GX I J , GX J I , G Y I J , G Y J I , GYY ) 

COMPLEX GI J,GNI J,GNJI ,GXI J, GXJ I , GY I J , GY J I , GYY 
DIMENSION COS AMU ( 200,25) , SI NAMU( 200,25) , AMU (200) , 
1C0EFG( 200) 

DIMENSION TEST1 (40 ) ,TEST2(40) , TEST3(40) , TEST30(40 ) , 

1 TEST4( 40 ) 

COMMON/ CON ST/H,0,DELTHE, SMI N,NPTS,NSPTS 

CCMMON/C PI / P I 

S UM 1 = 0 « 0 

SUM2=0 . 0 

SUM3=0 . 0 

SUM30=0 .0 

SUM4=0. 0 

DO 20 N = 1 » 40 

DO 15 K K = 1 ,5 

K = ( N- 1 ) *5 + KK 

SN I = S INAMU ( K , I ) 

SN J=S I NAMU ( K , J ) 

C S I =COS AMU ( K » I ) 

CSJ=CGSAMU( K, J ) 
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AK= AMU ( K) 

CCK=COEFG ( K ) 

VAL=AK*ABS(XV-XC ) 

IF (VAL.GE. 75.0) GO TO 10 
EXIJ=EXP(-VAL) 

GO TO 11 

10 EX I J = 0 . 0 

11 CONTINUE 

SUM1=SUM1+C0K*CSI*CSJ*EXI J 
SUM2=SUM2+CGK*CSI*CSJ*EXI J*AK 
SUM3=SUM3*C0K*AK*SNI*CS J*EX I J 
SUM30=SUM30+C0K*AK*SNJ*CSI*EX I J 
SUM4=SUM4+AK*AK*C0K*C SI *CSJ*EX I J 

15 CONTINUE 
TEST1 ( N ) = SUM1 
TEST2( N)=SUM2 
TEST3 ( N ) =SUM3 
TEST30(N)=SUM30 
TEST4! N ) =SUM4 

I F < N- 1 ) 20,20,16 

16 IF( ABS( SUM1) ,LE. 0.000001) GO TO 17 

I F ( ABS ( (TESTl(N)-TESTl(N-l) )/TESTl(N)).GT. 
1 TO 20 

17 I F(A8S( SUM30) .LE .0.000001 ) GO TO 18 

I F < ABS < !TEST30!N)-T5ST30!N-1) )/TEST30(N) ) . 
1 GO TO 2 0 

18 I F ( ABS ( SUM3) .LE. 0.000001) GO TO 19 
IF(ABS((TEST3(N)-TEST3(N-1))/TEST3(N)).GT. 

1T0 20 

19 IF(ABS! SUM2) .LE. 0.000001) GO TO 21 
IF(ABS( ( TEST2( N ) - TE ST2 ( N- 1 ) )/TEST2(N) ) .GT. 

1T0 20 

21 I F < ABS ( SUM4) .LE. 0.000001) GO TO 30 

IF! ABS! < TEST4! N)-TEST4(N-1 ) )/TEST4(N) ) .GT. 
1T0 20 
GO TO 30 

20 CONTINUE 

WR I TE ! 6 , 2 5 ) I , J 

25 FORMAT! 10X23 HG REE NS DID NOT CONV ERGE , 2X 2H I 
1 2X2HJ= I 2 ) 

30 CONTINUE 

IF (XV.LT.XC) S I GN=- 1.0 
IF (XV.GT.XC) S I GN=1 . 0 
IF (XV.EQ.XC) S I GN=0 . 0 

TERM4=4.*PI*CHYI*CHYJ*SIN!A*A6S< XV-XC) ) / < 2 
TERM5=4.*PI*CHYI *CHYJ*COS! A* ABS ( XV-XC) ) / ( 2 
TERM6=A*TERM5*SIGN 
TERM7=A*TE RM4* SIGN 

TERM8 = 4.*PI*A*Sr!YI*CHYJ*C0S! A*A3S( XV -XG ) ) / 
1SH2AH) 

TERM9=4.*PI*A*SHY I*CHY J*SIN ( A* ABS! XV-XC ) ) / 
1SH2AH) 

TERM80=4.*PI*A*SHYJ*CHYI*C0S( A* ABS! XV-XC ) ) 
1SH2AH) 

TERM90=4.*PI*A*SHYJ*CHYI*SIN(A*ABS(XV-XC ) ) 
1SH2AH) 

TERM10=A*A*TERM4 
TERM11=A*A*TERM5 
IF (J.E0.25) GO TO 70 
G I J=C MP LX { SUM1 + TE RM4,-T ERM5 ) 

GXI J= CMPLXl -SIGN* SUM 2+TERM6,TERM7) 

GY I J=CMPLX!-SUM3*TERM9,-TERM8) 

IF < I . EQ.25 ) GO TO 40 
GXJI=CMPLX(SIGN* SUM2-TERM6 , -TERM7 ) 
GYJI=CMPLX( -SUM30+TERM90, -TERM80) 

GNI J=CMPLX!-ANXI*SIGN*SUM2+ANX I*TERM6-ANYI 
1 TE RM9 ,ANXI*TERM7-ANYI*TERM8 ) 
GNJI=CMPLX!ANXJ*SI GN*SUM2-ANX J*T ERM6-ANY J* 
1TERM90,-ANXJ*TERM7-ANYJ*TERM80) 

GO TO 60 
40 CONTINUE 

GYY=CMPLX( -SUM4+TERM10 ,-TERMll ) 



.0010) GO 
GT. .0010 ) 
.0010) GO 
.0010) GO 
.0010) GO 

= 12 , 



. *A *H + S H2 AH) 
. *A*H + SH2AH) 

! 2 . *A*H +- 
( 2. *A*H + 

/ ( 2 .*A*H + 

/ ( 2 • *A * H + 



* SUM3*ANY I* 
SUM3G+ANY J* 
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60 CONTINUE 
GO TO 30 

70 G I J=CMPLX { -SUM! -TE RM4 , - TERM5 ) 

GNI J=CMPLX (ANXI*SIGN*SUM2-ANXI*TERM6+ANY I*SUM3-ANYI* 
LTERM9, ANXI*TERM7-ANYI *TERM8 ) 

80 CONTINUE 
RETURN 
END 



THIS SUBROUTINE INVERTS COMPLEX MATRICES TO SOLVE 
THE MATRIX EQUATION ALPHA! I , J ) *F ( I , 1 ) = HH<I,I) AND 
ALPHAU ,J)*F1(I,1) = PK (1*1$ 



SUBROUTINE COMAT ( N , M , A , B , 0 , I ) 

INTEGER C,H,R,Q,Z 

COMPLEX A, 6, 0, TT , P 

DIMENSION A(N,N) ,6(N,M) ,C( 100,3) 

D = ( L .0,0.0) 

DO 20 J = I , N 
20 C(J,3) = 0 

DO 21 K = I , N 
TT = (0.0, 0.0) 

T = 0.0 
DO 4 J = L,N 

IF (Cl J,3) .EQ. 1 ) GO TO 4 
DO 5 H = 1 ,N 
IF ( C ( H , 3 ) - 1) 15,5,12 

15 IF (T .GE. CABS( A( J,H) ) ) GO TO 5 

R = J 
Q = H 

T = CABS ( A ( J , H) ) 

5 CONTINUE 
4 CONTINUE 

C (Q,3) « C( Q,3) ♦ 1 
C ( K , 1 ) = R 
C ( K , 2 ) = Q 

IF (R . EQ. Q) GO TO 11 
D = -D 

DO 8 L = 1 ,N 
TT = A ( R , L ) 

A ( R ,L ) = A ( Q , L ) 

8 A(Q,L) = TT 

I F ( M .LE. 0) 30 TO 11 
DO 2 L = 1 ,M 
TT = B ( R , L 1 
B ( R »L ) = B ( Q , L ) 

2 B(Q,L) = TT 
11 P = A(Q,Q) 

A ( Q ,Q) = (1.0, 0.0) 

DO 13 L = 1 , N 
13 A ( Q , L ) = A(Q,L)/P 

IF (M .LE. 0) GO TO 1 
00 3 L = 1,M 

3 B ( Q , L J = B ( Q , L ) / P 
1 DO 21 Z = 1 , N 

IF (Z . EQ. Q) GO TO 21 
TT = A ( Z , Q ) 

A(Z,Q) = (0. 0,0.0) 



DO 16 L = 1 , N 

16 A ( Z , L ) = A ( Z , L ) - A ( Q . L ) *TT 
IF (M . LE. 0) GO TO 21 

DO 17 L = 1 , M 

17 B ( Z , L ) = B ( Z , L ) - B(Q,L)*TT 
21 CONTINUE 

DO 19 II = 1 , N 
L = N t 1 - II 

I F { C ( L , 1) .EQ. C ( L , 2 ) ) GO TO 19 
R = C ( L , 1 ) 

Q = Cl L » 2 ) 
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DO 7 K = 1,N 
TT = A ( K »R ) 

A ( K » R) = A(K,Q) 

7 A l K »Q ) = TT 
19 CONTINUE 

DO 18 K = 1,N 

IF ( C ( K » 3 ) .NE. 1) GO TO 12 
18 CONTINUE 
I = 1 
50 RETURN 
12 I = 2 
GO TO 50 
END 



/ /GO. $ Y SI N DO * 
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